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Chapter 1 

Introduction 



1 .1 Overview and Notation 

We begin with an overview of filters, filter banks, and wavelets. We want to indicate, first in 
rough outline and then in detail, the connections between these three topics. Our immediate pur- 
pose is to open up the problem and the language— starting with the filter coefficients A(n). The 
choice of those coefficients is the crucial decision. Their properties govern all that follows. 
Each step is a natural development from the one before: 

(1) A filter is a linear time-invariant operator. It acts on input vectors x. The output vec- 
tor y is the convolution of x with a fixed vector ft. The vector h contains the filter coefficients 
h(0), A(l), ft(2), .... Our filters are digital, not analog, so the coefficients h(n) come at discrete 
times t = nT. The sampling period T is assumed to be 1 here. The inputs x(n) and outputs y(n) 
come at all times t = 0, ±1, ±2, . . .: 

y(n) = ^2 h{k)x(n -k) = convolution A * x in the time domain. 

k 

One input x = (. . . , 0, 1, 0, . . .) has special importance — a unit impulse at time zero. The 
input has x(n — k) = 0 except when n = k. The sum in the convolution has only one term, and 
that term is h(n). This output.? (n) = h(n) is the response at time n to the unit impulse x(0) = 1, 
It is the impulse response A(0), A(l), , . . . , h(N). 

In a moment the same filter will be described in the frequency domain. Convolution with the 
vector h will become multiplication by a function H. It is the simplicity of multiplication that 
makes this subject a success. The action of a filter in time and frequency is the foundation on 
which signal processing is built. 

(2) A filter bank is a set of filters. The analysis bank often has two filters, lowpass and high- 
pass. They separate the input signal into frequency bands. Those subsignals can be compressed 
much more efficiently than the original signal. Then they can be transmitted or stored. We are 
describing "subband coding" and its applications. At any time the signals can be recombined 
(by the synthesis bank). 

i ;„■ ft is not necessary to preserve the full outputs from the analysis filters. Normally they are 
^ cbwnsampled. We keep only the even components of the lowpass and highpass filter outputs. 
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1.1 Overview and Notation 



expansion displayed in equation (1.1) by w JK (t) and integrate: 



I 

J-o 



f(t)w ]K {t)dt 



=b JK r 



(w JK (t)) 2 dt. 



(1.3) 



All other terms in the sum disappear because of orthogonality. Equation (1.2) eliminates all in- 
tegrals of Wjt times wjk , except the one term that has j = J and k = K. That term produces 
(u)jk(?)) 2 - b JK is the ratio of the two integrals in equation (1.3). 

As we describe the connection between filter banks and wavelets, you will see that it is the 
"highpassfilter" that leads to 10(f). Hie "lowpass filtef leads to a scaling function </> (f). In most 
constructions the lowpass filter comes first — the scaling function is obtained before the wave- 
let. In fact the scaling function (in continuous time) comes from infinite repetition LL ... L of 
the lowpass filter, with rescaling at each iteration. The wavelet follows from <f>(t) by just one 
application of the highpass filter. 

Multiresolution 

At a given resolution of a signal or an image, the scaling functions <f> (2 J t - k) are a basis for the 
set of signals. The level is set by j 9 and the time steps at that level are 2~K The new details at 
level j are represented by the wavelets w (2 J t — fc). Then the smooth signal plus the details, the 

t 0's plus the iu's, combine into a multiresolution of the signal at the finer level ; + 1. Averages 

|) come from the scaling functions, details come from the wavelets: 

signal at level j (local averages) 

+ signal at level y H- 1 

details at level j (local differences) f 

; That is multiresolution for one signal. When we apply it to all signals, we have multiresolution 
for spaces of functions: 



Vj = scaling space at level j \^ 



MM 



Wj = wavelet space at level j f 



V)+i = scaling space at level j + 1 



p This idea of multiresolution is absolutely basic to wavelet analysis. Again, we are only intro- 
• during it We are sending a coarse signal to the reader, not the details. You only have the input 
gatjevel 1. 

Thus the signal is divided into different scales of resolution, rather than different frequen- 
|: h cies * "time-scale plane" takes the place for wavelets that the "time-frequency plane" takes 
fillf ft* filters. Multiresolution divides the frequencies into octave bands, from co to 2a>, instead of 
JhB^a I «^fonn bands from co to co + Aco. The compression of a graph, when f(t) is replaced by f (2t), 
I ; ' S ex Pansion of its Fourier transform from F(co) to \ F (f ). Frequencies shift upward by an 

'^^l^octe^^hen time is rescaled by two. You will see how the time-frequency plane is partitioned 

naturally into rectangles of constant area (Figure 1.1). 
yv e y/ x : matching of long time with low frequency and short time with high frequency occurs 
^^^^;2 1 natural way for wavelets. It is one of the attractions of a wavelet decomposition. 
Iift&^v%^J^ er: We reprinted in Appendix A an article on wavelets published in the 
^iiMi^^rf^ Scientist of May 1994. This article introduces wavelet notation through its correspon- 
^^P^^^P^* musica l notation. In music, each note specifies a frequency and a position in time. 
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Multiresolution will be described first for subspaces Vj and Wj. The ^ 
• apt, ^eiet snace W, is the difference between V; and V J+X . The sum oj v, 

increasing. The wavelet space w s is. jj translation to t - k define 

Wj is V i+1 . Then these extra conditions involving ddation to It and translation 

* If /ft) is in V, then /ft) and /(2r) and all /ft - *) and /(2r - k) are in ^i- 

In theend, one wavelet generates a whole basis. ^^.^-^^iS.*-: 
translation (all j and all k). There are six steps toward this goal, and we take them 

1. An increasing sequence of subspaces Vj (complete in L 2 ) 
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2. The wavelet subspace Wj that gives Vy + Wj = V i+ i 

3. The dilation requirement from /(f) in V) to f(2t) in V /+ i 

4. The basis <£(r - *) for y o and w(t - k) for W 0 

5. The basis 0(2' f - *) for Vj and w(2 ; f - k) for W y 

6. The basis of all wavelets w(2 j t - k) for the whole space L 2 . 

shortcut to multiresolution. Before those six steps, may I mention one shortcut step that 
trts with the filter coefficients h(k). That step is to solve the dilation equation for the scaling 
notion <t>(t): jgmmMBSK < ; ^ 

|*(f) = E2M*)*(2r-*)-| 

ie first requirements on the coefficients are £fc(fc) = 1 and £ (-1)* A(4) - 0. The full 
quirement is Condition E in Section 7.2. When this is satisfied, <f>(t) can be computed. Then 
H2 J t - k) } is a basis for V,- . These spaces are automatically increasing and complete and shift- 
variant and connected by dilation, thus multiresolution is achieved. 

Scale of Subspaces 

ach Vj is contained in the next subspace V j+X . A function in one subspace is in all the higher 
iner) subspaces: 

VoCViC-'-CVjCVj+iC--- 
. function /(f) in the whole space has a piece in each subspace. Those pieces contain more and 
tore of the fall information in /(f). The piece in Vj is />(f). One requirement on the sequence 
f subspaces is completeness: 

fj(t) -> /(f) as j oo. 

The first example will not have the dilation feature required for multiresolution: 

Sample 6.1. V) contains all trigonometric polynomials of degree < /. 
lertainly Vj is contained in V j+X . The spaces are growing. (Since Daubechies uses - ; where we 
se /, her subspaces are decreasing. Most authors now use an increasing sequence, for simpler 
lumbering.) The piece of f(t) in Vj is the partial sum fj (f ) of its Fourier series: 

fj(t) = c * eih is * e P iece in V l" 

rhis is the projection of /(f) onto Vy. The exponentials e ikt are orthogonal, so the energy in 
fjO) is the sum of \c k \ 2 over low frequencies \k\ < j. The energy in /(f) - fj(t) is the sum 
>ver high frequencies \k\ > j. This approaches zero as ; -> oo. Therefore the sequence Vj is 
complete in the whole 27r -periodic space L 2 . 

Now we identify the second family of subspaces. Wj contains the new information 
&fj(t) = /j+i(f) - //(f). This is the "detail" at level From the viewpoint of individual 
functions, 

|/)(0 + A//(0 = /)+i(0- } (6-D 
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From the viewpoint of the subspaces they lie in, this is 

Vy0Wy = Vy +1 . 

Each function in Vy+i is the sum of two orthogonal parts, fj in Vy and A/y in Wy . In our e? 
pie, the new information A/y = f j+ \ - fj is the new term that enters /y + i (t): 

A/ J (0 = c i+ie '0 + ». + c _._ ie -.a + .)*. 

The space Wj contains terms of exact degree j + l. Those are orthogonal to all terms of de 
< j. Together, these orthogonal complements Wy and Vy produce Vj+\. This is an impo 
part of multiresolution. 

The spaces Wj are differences between the Vy. 
The spaces Vy are sums of the Wj. 

We can call Vy a partial sum by recognizing how the W's add to V. Start from 

Vb©Wb«V, and Ki0JVi = V2. 

Substituting the first into the second, V 2 is the sum of three mutually orthogonal subspaces 

Vo0W o 0W, = V 2 . 

When you add details up to and including Wy , you have Vy +l : 

Vb0Wo0W!0..-0Wy=Vy +l . 

For the functions in those subspaces, this equation is just 

/o(0 + A/ 0 (0 + A/,(0 + - - • + A/y(0 = /y + i(f). 

That sum is / 0 + (fx - /o) + Cfe - f\) + • • • + (fj+x - fj)- It 'telescopes" into /y +1 . 

In practice, we can construct the spaces Vy and take differences. Or we can construe 
Wj and take sums. It is like differentiating an integral or integrating a derivative. Starting : 
a basis is like starting from the Wj\ then the Vy are partial sums as above. For the expone 
basis e ik \ this works fine. It also works for Haar wavelets, which are so simple. But for 
wavelets, that direct search for w(t) was very difficult. 

Important The construction of wavelets has succeeded by finding the Vy first. We begin 
the scaling function 0(f), not the wavelet! Its translates <f>(t — k) go into Vo. Rescaling t< 
gives Vy. Then the wavelet spaces Wj are the differences between Vy+i and Vy. The fiinc 
in Wj are the details at the yth scale. 

Similarly for filter banks, we design the lowpass filter by choosing c(k). Then the high 
coefficients d(k) are easy. To maintain this analogy between continuous and discrete, we dr 
"logarithmic tree" of filter banks. At each step, the highpass filter (with downsampling as u: 
produces the detail A/y in Wy. The space Vy+i is the sum of Vo, W 0 , . . ., Wy. 

It is convenient if Wj is orthogonal to Vy. Each Wj is then automatically orthogonal t 
other Wjt. Reason: Uk < j, then W* is contained in Vy — which is perpendicular to Wy. 
completeness condition can be restated as 

7=0 
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With orthogonaUty of each piece /,(,) to the next detail A/y (,), these subspaces are orthogonal 
But we emphasize now that orthogonality is not essential. eortnogonal. 

A nonorthogonal example comes direcdy from any nonorthogonal basis b 0 (t), b,(t) 
The piece /;(;) includes all the terms through bj(t): W 

Sum up toy: /y(«) « dA(f) «in V. 
Next term: A/y(f) = c j+l b j+l (t) is in IVy. 

^L atem K i ?° tlOSt ' jUSttheOrth0g0nality - DCW *** V >« isstillthe-ATm^-of V- 
and Wj which intersect only at the zero vector. The angle between subspaces can be less than 
90°, as long as every f j+l in V m has exacdy one splitting into fj + Afj : 

;: Vyn^ = {0 } and V j + W j = V j+1 \ (6fi) 

This nonorthogonal situation applies to biortkogonal filters and wavelets (Section 6 5^ TW 
Wj is orthogonal to a different subspace Vj . The extra freedom can be put to^ use 







v j + l 







— v j 


v j + l 





Figure 6.1: An orthogonal sum and a direct sum. Both written VjQWj = Vj+l and both aUowed. 



The Dilation Requirement 

So far we have an increasing and complete scale of spaces. Each V, is contained in the next 
Vy+i. For multrresolutton, the crucial word scale carries an additional meaning. V M consists 
of all rescaled functions in V } : -S- "y+i consists 

Dilation: /(f) is in V,- /(2f) is in V, + , . 

25 SSSfj-j ? } ! hang « T ^ 35 ^ 38 ^ ° f /W - ° n a ma P' * e «* " doubled 
that '^^'^^ ' 16 ° S a ^ a ^ e " ^ 6,000,000 : 1 its height is a half page. The length 

that represents a mile is cut in half. This length is At or Ax or A "»length 

The example using /(,) = c Jje~>J' + . . . + ^ does not meet ^ rescali requirement 

Se^^ 

toco^ nT nCy e T meS2y - Thenew space V0 + , ^required 

wteveTari ^. neW T CnC,eS - T ° SatiSfy ^ the partial sumTgo an 



Multiresolution example : 
Next detail : 



fj(*) = £e t «* for |*| <2' 
A /K0 = T.c k e ikt for 2> < |*| < 2^' +I . 
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This is a genuine multiresoiution, in which Vj and Wj have roughly the same dimension. It is 
the Litdewood-Paley decomposition of a Fourier series* into octaves instead of single terms. 

This is a chief part of the mathematical background. To fit the requirements precisely, when 
/(f) is defined on the whole line -oo < f < oo, we should use all frequencies co and not just 
integers: 

tot J\<o\<V 

Now the spaces V) go down the scale toward j = -oo, as well as up the scale. The continuous 
frequency co can be halved as well as doubled. The basis functions become sine functions, by 
the sampling theorem. Continuous frequency but discrete basis, as is normal for L 2 . And the 
nested spaces include j < 0: 



• • • C V-i C V 0 C • ■ • C Vj C V>+i C • • • 
In addition to completeness as j -> oo, we require emptiness as j — oo: 



P| Vj = {0} and (J Vj = whole space. 



(6.7) 



(6.8) 



Emptiness means that J //(f) II 0 as 7 -> -oo. Completeness still means that fj(t) /(f) 
as j oo. The detail Afj == fj+\ — fj belongs to Wj and we still have 



(6.9) 



This can be an orthogonal sum, with A/) orthogonal to fj. It must be a direct sum, with 
Vj f) Wj = {0}. The reconstruction of /(f) from its details Afj can start at j = 0 as before, or 
it can start at 7 = — oo: 

/(0 = /o(0 + f>./)(f) or /(*) = £ A// (0- 

0 -oo 

The sum of subspaces can start at 7 = 0 or j = -00. When the sum stops at / > 0, we have 
the subspace V)+i: 

y j 
V J+l = V Q + J2Wj or V J+l = 2 W,.. 

The left sum includes the scaling functions in V 0 ._The sum on the right involves only the wave- 
lets. That form includes all the very large time scales Af = 2~ j as j -> -00. 

In practice we use the first sum. Our calculations begin at some unit scale. The scaling func- 
tions at j = 0 and the wavelets with j > 0 are the basis. I suppose the scaling functions at level 
7 = 7 and the wavelets with j > J are another basis. 



The Translation Requirement and the Basis 

Instead of rescaling /(f), we now shift its graph. This is translation, and it leads to the funda- 
mental requirement of time-invariance in signal processing. The subspaces are shift-invariant. 

If fj(?) is in Vj then so are all its panslates fj(t — k). 
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Suppose f(t) is in V 0 . Then f(2t) is in and so is ft>t m o • ^ 

With translation we are committed to working on the whole line ,sn ^ , 
riodicity. A particular /(,) may have compact suppo£ * £ Xle^e V^'fT ^ 
together) is shift-invariant For finite intervals the ^2 f ! ? T ° 1 fanCt, ° nS 

the whole space V Q . They may even be ormono^,. ^ translates ^ S P*> 

later, is that V 0 contains suZ^n ^ """^ to ta Weakened 

^reexistsmsothatm-k^isanorthonormalbasisforVo. 

When the functions ^(r - t) are an orthonormal basis for V n the re <™w fi ^ 
^(^-^wiUbeanorthonormalbasisforV, AtscalinelZl, A \ teSCaled 
| -nonnalizedbyZ^.WecoUectaUme Jui^lS 

| 

| Multiresolution Analysis 

j| 7fc subspaces V, iorfsfr requirements 1 to 4: 

I 1. ^/CV y+1 and nVy = {0) and U^ = L 2 (completeness). 

2. Scafe invariance: f(t) € V, /(2/) <= 7;+,. 
|- 3. Shift invariance: f(t)eV 0 *=*f(t-k)e V 0 . 
■ 4> Shi fi-invariant basis: V 0 has an orthonormal basis {0( f - A:)}. 
-V 4 - ^-«vor^W: V 0 has a stable basis (Riesz basis) {<t>(f-k)}. 

1 i CO 

iP 7 > w - fl v*0y* (0 is the piece in V; . 
h the orthogonal case, the energy in this piece is 



: ' ■ CO 



IP: lly)l|2= i S 0 1 ^ 12 - (6.10) 

H| ^ at least three ways to construct or describe a multiresolution- 
° :f;v * y * e spaces V) 
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2. By the scaling function 0(0 

3. By the coefficients 2fc(Jfc) in the dilation equation. 

i/ ««^tHAirhaQP<; Then we move to description 3 and the 



come fet THen we study the mMon amnion, ami const.** wavelets. 

_^ A. 



fWinVj m-DinVj f(2t)inV j+1 

Figure 62: TMdoo stays in V,. Dilation moves into V, + „ Why is /(») - /(0 in Wfl 



Examples of Multiresolution 

/(f) = /(integer part of f). 
Ue mnrfon ,<2» in V, is then conaan, on ^^ZSZ 

wise constant The step from j to j + 1 rescales time by 2 and produces V j+1 . What about 
basis? The simplest choice is the box function: 

J 1 for0<f<l is ort h 0 gonal to its translates 0(t-*)- 
0 otherwise 

Every function in V„ is a combination of boxes /(*) = £/(»)*(>-«>• So requirement 4 i 
satisfied by the box function 0 (t). 

^ rwinnnns oiecewise linear functions. The functions /(.) are now linear between eac 
The spaces are right for multiresolution. Is there a shift-mvariant basis? 
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The basis function that comes to mind is the hat Junction H(t), equal to one at t — 1, 
and linear between its values H(n) = 6(n - 1). The translates H(t - k) generate all piece- 
wise linear functions on unit intervals. Any function /(/) in V 0 can be expressed as 
Y, f(n - l)H(t - n). However H(t) is not orthogonal to the neighboring hat H(t - 1). The 
product H(t)H(t - 1) is positive on the one interval 1 < * < 2 where the hats overlap. Its 
integral (the inner product of the hats) is not zero. 

We must work harder to find an orthogonal basis, and the eventual <f> (t) will not have compact 
support. Or else we keep this non-orthogonal basis. 

3. Discontinuous piecewise linear functions. Now / (t ) in Vb may have a jump at each mesh- 
point t = n. There is a value /(*_) from the left and a value /(n + ) from the right. The hat func- 
tion is still in the space, but so is the box function! The spaces Vj are clearly shift-invariant and 
scale-invariant. If f(t) is linear between integers (where it jumps), then f(2t) is linear between 
half-integers (where it jumps). 

There are two degrees of freedom at each meshpoint, the values f(nJ) and /(n+). Therefore 
two scaling functions <f>\(t) and <fe(r) are required for a shift-invariant basis. They can both be 
supported on the unit interval, and they can be orthogonal: 

0i (0 = box function and <k(t) = sloping line = 1 - 2t. 

The union of {<f>\(t — k)} and {faQ - k)} is an orthonormal basis — which illustrates the idea 
behind "multiwavelets". The usual dilation equation for 4> (0 becomes a vector equation for tp\ (f ) 
and <h(0- The coefficients c(k) in that equation are 2 x 2 matrices. The associated filter bank 
in Section 7.5 contains "multifilters". 

4. Cubic splines. Vq consists of piecewise cubic polynomials on unit intervals, with f(t ) and 
f'(t) and /"(0 continuous. The third derivative /"'(f) may jump at the integers t = n, so 
the cubics are different in neighboring intervals. We have shift-invariance and scale-invariance, 
when Vi contains the cubic splines on half-intervals. This is the main point: Approximating 
subspaces on regular meshes automatically fit the requirements for multiresolution. 

The shortest cubic spline is a B-spline. It consists of different third-degree polynomials on 
the four unit intervals within 0 < t < 4. The letter B stands for basis, but not for orthogonal ba- 
sis. In complete analogy with the hat function, which is a linear spline, the scaling function <f>(t) 
for the cubic splines cannot have compact support if we insist on orthogonality. An orthogonal 
basis [<f> (t — k) } does exist, but it requires Fourier analysis to find it. 

The cubic B-spline satisfies a dilation equation with very simple coefficients, proportional to 
1,4,6,4,1. But those coefficients do not lead to orthogonal filters. We can stay with these co- 
efficients and go to biorthogonal filters (the best plan). Or we can orthogonalize, losing compact 
support and reaching a filter with infinitely many coefficients. Section 7.4 develops the theory 
of splines. 

5. Daubechies functions. The search for orthogonal filter banks leads to the four coefficients 
of a "maxflat" lowpass filter. The response C(co) has a double zero at the highest frequency 
<*> = Jr. This is maximal flatness, with four coefficients and orthogonality: 



c(0), c(l), c(2), c(3) = 1 + V3, 3 + >/3, 3 - >/3, 1 - times 1/4^2. 
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Their sum is \/2. Their sum of squares is unity. Uiey are orthogonal to their double shifts, be- 
cause c(0)c(2) +c(l)c(3) = 0. From these coefficients Daubechies constructed <f>(t) by solving 
the dilation equation 

_ 3 

0(O = A/2£c(*)*(2f-*). 

The solution comes in the next section. The zeroth space V 0 contains every <f>{t - k). Those 
functions are an orthononnal basis. The rescaled functions <f>(2 J t — k) span Vj. 

This is our best description of the Daubechies spaces V}, to give the dilation equation for 
In Examples 1-4, we started with the spaces. In Example 5, Daubechies started with 
the coefficients and found <t>{t)— which produces the spaces. Either way, we have the scale- 
invariance and shift-invariance of muitiresolution analysis. 

The Dilation Equation 

The space Vb is contained in V\. Therefore <f>(t) is also in Vi. It must be a combination of the 
basis functions 2 X/ 2 <p (2t—k) for that subspace. The coefficients in the combination will be called 
c(k). Bring the factor 2 1/2 = outside: 

1 VbC V! means <f>(t) = c(k) 0(2f - k). 1 (6 U) 

This is the dilation equation. It is a two-scale equation, involving t and It. It is also called a 
refinement equation, because it displays <p(t) in the refined space Vi. That space has the finer 
scale At = 1/2, and it contains <f>(t) which has scale At = 1. 

To emphasize: The dilation equation is a direct consequence of VqCV\. It is not an extra 
requirement! There will be a finite set of coefficients c(0), . . . , c(N) when <f>(t) is supported on 
[0, N]. In general, <f>(t) has infinite support and we need infinitely many c(n). 

Tofindc(n), multiply the dilation equation (6. 11) by J2$(2t-n). Integrate and use orthog- 
onality: 

>/2 f* WWQt-n)dt~c(fl). (6.12) 

If <p(t) is the unit box and <l>(2t) is the half-box, this gives c(0) = V2/2 and c(l) = V2/2. The 
dilation equation for the box function then has coefficients 1 and 1: 

0(O = <£(2O + *(2f-l). ( 6 - 13 > 

From orthogonality of the basis {<p(t —A:)} we have double-shift orthogonality of the dilation 
coefficients c(Jfc). And unit energy in 4>{t) gives a unit vector of c's: 

Double-shift: ^c(k)c(k - 2m) = 5(m). Unit vector: £ |<r(*)| 2 = 1 (6.14) 

For proof, multiply the dilation equations for 0(f) and <f>(t — m) and integrate. Orthononnality 
of the yields double-shift orthogonality of the c's: 

f°° <t>(t)<t>(t - m)<fr = £c(*)c(* - 2m ) = 5 < m )- < 615 ) 

J-oo * 
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The coefficients c(k) go into an orthonormal filter bank! Starting with the spaces V) in a multi- 
resolution, the dilation equation has brought us back to filters — where the key matrix is 
L = (4 2)C. Double-shift orthogonality becomes LL T — /. The rows of L contain the dou- 
ble shifts Lij = c(2i - ;). 



Box example: 



J2 



1 1 



1 1 



1 1 



Daubechies example: 



L = 



c(3) c(2) c(l)< c(0) 

c(3) c(2) c(0) 
c(3) c(2) 



To end this section, we have to identify the wavelet spaces Wj. 



The Wavelet Equation 

The scaling functions # (2 s t - k) are orthogonal at each scale separately. But <p (t) is not orthogo- 
nal to <f>(2t). They are not orthogonal across scales; the level j must be fixed. The function 0(0 
in Vo is also in V\ (the dilation equation). Orthogonality across scales comes firom the wavelet 
subspaces Wj and their basis functions w jk (t). We study those now, from three starting-points: 

1. The spaces W). 2. The wavelets w(t). 3. The coefficients d(k). 

Use Method 1 if you have the Vj . Their differences yield the spaces Wj . Use Method 2 if you 
can identify the wavelets. Just shift and rescale. Use Method 3 if you have the numbers c(k). 
The alternating flip yields = (— l)*c(N — fc). Then w(t) comes from the wavelet equation 
below, and Wj contains the combinations of w(2 j t — k). 

The box function gives an example in which all three approaches will work. We construct 
W 0 and w(t) and the <f s: 

1. From the subspaces: V 0 contains constant functions on unit intervals, and V\ contains 
constant functions on /w//-intervals. The space Wq is in Vi (therefore constant on half- 
intervals). It is orthogonal to Vb, so the integral over each full interval is zero. This fact pro- 
duces the complementary subspace Wb, orthogonal to Vb inside V\\ 

W 0 = { constants on half-intervals with /(n) + f{n + 1/2) = 0 }. 

Vb © Wo does give Vi . Combining equal values at n and n + 1 /2 from Vb with opposite values 
from Wq gives any values f(n) and f(n + 1/2) for V\. 

2. From the wavelets: The important function in Wq is the up-down square wave: 



Haar wavelet w (f) = 



1 forO<f<| 
-1 for \ < t < 1 
0 otherwise. 
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0 ± 



1 ± 



Figure 63: Two bases for V,: Halfsize boxes 0(2r - k) or Ml boxes 0(r - *) plus uptown 
Haar wavelets w(t-k). 

This is orthogonal to the box function 0(0- It is orthogonal to translates of 0 and also to its own 
translates (there is no overlap of »(f ) with w(t - 1)). More than that, multiresolution says that 
the wavelet w(f) is orthogonal to rescalings of itself and to translates ofrescalmgs: 

fCO 



f°° w(t)w(2 j t -k)dt = 0 unless j = k = 0. 



The translates of u>(0 span Hfc. H- translates of w&t) span W) Those wavek t«-»»* 
thogonalbecauselVo C V, and V, X W,. (Exchange; and 0 if j is negate.) Fro. o^ogonj 
spaces we have orthogonal basis functions. Then completeness makes the whole orthonormal 
system {2 J/2 w(2h — k)} a basis for L 2 . 

3. From the coefficients c(0) = c(l) = l/T* The flip construction gives d(0) = 1/72 and 
dm - -1/72 Those coefficients go into the wavelet equation: 

*-*).( (6-16) 

This equation pr^^Se^e^SnTS^liS functions » 
solveX ThVwavelet is u>(0 = 0(2/) - 0(2f - D- This is a half-box minus a shifted half-box. It 
is the uo-down square wave, which is Haar's wavelet 

OuSexaVf^ 
the four d's (to normalize, divide again by 472): 

<f(0), rf(l), d(2), d(3) = 1 - 73, -(3 - 73), 3 + 75, -(1 + 73). 

Tbeir sum is zero. Their sum of squares (normalized) is 1 . They axe orthogonal to .Ad. : double 
shifts because the c's are. The wavelet equation gives the Daubechies wavelet rvjt), whic xtm 
noSip'fonnula. The orthogonality to u,(r -*) and<Ht-k) is only known Meetly ~-f^ 
redTle-shiftorthogonalityofthecfs. Tne stnicture of multiresolution gives crucial mforma- 
Hon that we cannot find in a table of integrals. ^; Q ^i v 
The actual construction of 0(0 and u,(r), and the drawing of their graphs, is immediately 

ahead. 

Examples (Strange but beautiful.) Suppose 0(0 is the ^elta ^f 0 ^]^^^ 
butcontinue anyway. The space V 0 contains combinations E<z(n)S(r -n) of Mut functions at 
the integers. What orthogonal wavelet w(t) goes with this scaling function 5(0? 

By scale invariance, V, contains5(2r-n). The spikes for V, areaU = 0, ± 5 , ±1, • - , 
V 0 holds the delta functions at integers, W 0 contains delta functions at the midpoints t-n + r 
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lie integers and the midpoints combine to give V 0 © W 0 = V\ . The wavelet is the delta function 
" = 2- 

Similarly, V) contains delta functions at t = n/2 j . Wj contains delta functions at the mid- 
oints (n + \)/2>. What is W- X ? Its delta functions are at t = (n + \)/2" x = 2n + 1. These 

re odd integers ±1, ±3, ±5 The spacing between them is 2 as expected. Then WL 2 has 

elta functions at ±2, ±6, ±10, ... with spacing 4. The union of all Wj has delta functions at 
U binary points. 

The dilation equation for the delta function is S(t) = 28(2t). The only nonzero coefficient 
; h(0) = 1. The filter is the identity. The wavelet equation with only one term is w(t) = 
8(2t - 1) = S(t - \). This confirms what we found, that W 0 contains delta functions at all 
lidpoints between integers. Notice! An odd number of coefficients (one) means that N = 0 
sven). The alternating flip must shift by an odd integer, for double-shift orthogonality. So the 
onzero highpass coefficient was d(l) notrf(O). 

To linger one last second on this trivial great example, the double-shift matrices from the low 
nd high channels are L = (|2) and B = ( j 2)( delay ). 



x(n) 



-1 



(2 



-1 



r I 2 



L|2 



(x(0),x(4),x(8) ) 

(x(-2),x(2),x(6) ) 

(x(-1),x(1),x(3), ) 



Figure 6.4: The lazy filter H = / leads to delta functions. 



lie Scaling Function is Supported on [0, N] 

l remarkable feature of 4>(t) is that it is zero outside the interval 0<t<N. This could never 
|ppen to a one-scale difference or differential equation (homogeneous). The solutions would 
| combinations of X n and e u 9 and only occasionally zero. The compact support of 4>{t) comes 
nom the two scales in the dilation equation 



N 

^(o = £>(*) 0(2/-*). 

*=0 



(6.17) 



lieorem 6.1 The scaling function <f>(t) is supported on the interval [0, N], 



^of. Suppose we know that the support is aflnite interval [a, *]. Then <f>(2t) is supported on 
2»-2l- The shifted function <p(2t - k) is supported on [*±£ t The index k goes from zero 
? jy , so the right side of the dilation equation is supported between § and Comparing with 

leftside, - 

[a b + Nl 
2 leads to a = 0 and b = N. 



6.2 Wavelets from Filters 



187 



This section will do part of each job, the recursive part. This shows how multiresolution (for 
functions) connects to subband filtering (for vectors). The three parts that we can do immediately 
are: 

1. Compute ajk and bjt recursively from a;+i f * (and vice versa). 

2. Set up a recursion (the cascade algorithm) to construct <p(t). 

3. Prove orthogonality for 4>jt(f) and wjk(t) from orthogonality of c's and <f s. 

Those are the three subsections. Later we have to initialize the recursion in 1, execute and study 
the cascade algorithm in 2, and derive other properties in 3. 

Wavelet Coefficients by Recursion 

Suppose f\ (t) is in V\ . It is a combination of the basis functions ^/2</>(2t — k). These functions 
<t>\k(t) are at level 1. Multiresolution splits this level into Vi = Vb © W 0 , so f x (f) is also a 
combination of the basis functions for Vb and Wb. Those basis functions are <f>ok(t) = 4>(t — k) 
and u)Qk(t) = w(t — k): 

E «u 0i*(O = E *Q* fad) + E *ot WQ*(0 

(6.19) 

= E*t* tf('-*0 + E*o* 

We are computing a change of basis. Given the coefficients a\k(t) in the V\ basis, we want the 
coefficients ao* and fro* in the Vo © W 0 basis. The same step will apply at every level. It takes 
us from the coefficients a/+i,jt in the basis for Vj+u to the coefficients and bjt in the bases 
for Vj and Wy. This is the recursion that makes the wavelet transform fast. 

We will suppose that these bases are orihonormal. Later in this section we prove this prop- 
erty (assuming the cascade algorithm uses orthogonal filters and converges). Orthonormality 
makes the formulas easy and it makes the inverse easy. Section 6.5 will derive the biorthogonal 
recursion, when orthogonality is not assumed. 

To find the recursion, shift equation (6. 1 8) by k and set n = £ — 2k: 

Dilation equation: <f>(t - k) = E a/2c(/i) 0(2* - 2k - n) = E*tf - W *ii (0 

Multiply by f\ (t) and integrate with respect to t . Since the basis functions are orthonormal, the 
integral gives the coefficients of f\ (t) in each basis: 

a ok =Y t c(i-2k)a H and = ^d(t -2k)a xt . (6.21) 

. . This is the key recursion. It is the action of a filter bank, which inputs an and outputs oq* and 

# &o*. But we have to watch indices, because an ordinary convolution would be E c (* — t) &u 

|;; and downsampling would give E c (2* ™ ^) fl u- There is a time reversal between this filter C 

| and the transpose filter C r that appears in the recursion (6.21): 

c r (n) =<?(-/») and d T (n)=d(-n). (6.22) 
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How do we know that the support is a finite interval in the first place? From the cascade 
algorithm. The box function 0 (O) (O is supported on [0, 1], When this box is substituted into 
the right side of the dilation equation, the function <£ {1) (r) that comes out has support [0, 
Then <f> 0) (t) is substituted into the right side and the result <£ (2) (f ) is zero outside [0, . The 
limiting function #0) is certain to be zero outside [0, N], This cascade is studied in the next 
section. 

It will be useful to reach the same conclusion based on the Fourier transform (Section 6.4). 
That argument can assume less about the filter coefficients. We mention that there are never gaps 
where 0(0 is zero on an interval inside [0, Afl. And if the highpass coefficients h x (k) run from 
k = 0 to k = N, then the wavelet w(t) = £ 2k x (k)4>(2t - k) has support [0, \(N + AT)]. The 
last term <j>(2t - N) is zero after 2t - N reaches N. 

Problem Set 6.1 

1. Explain why the scaling requirement, that ^(f) is in Vjjf and onl^if /(2f) is in V/ + i, can 
be restated as /(a>) »s in Vj if and only if f(2o)) is in V,_i. Here Vj is the space of Fourier 
transforms of functions in Vj. 

2. For the space Vb of piecewise constant functions in Example 1, show that the only shift- 
invariant basis (p(t - k) contains box functions. What is the corresponding statement about 
allpass FIR filters? 

3. For piecewise constants, show that /(f) is in L 2 if and only if f(n) is in I 2 . 

4. Find 2 by 2 matrices c(0)andc(l) so that the box function <fo(f) and sloping linefoO) = l-2f 
in Example 3 satisfy 

[tSMiSM J£S=S ]- 

5. If /(f) is in V Q and g(f) is in V u why is it generally false that g(f) - f(t) is in W x ? 

6. What multiresolution requirements are violated if Wj consists of all multiples of cos(2'f)? 



6-2 Wavelets from Filters 

The previous section reached the dilation equation and the wavelet equation: 

<£(f) = £V2c(n)0(2f -n) and u>(f) = £V2 d(n) <t>(2t - n). (6.18) 

Those equations are the crucial connections between wavelets and filters. Historically, their de- 
velopment was separate. Now you have to see them together. The lowpass filter c(0) c(N) 

determines the scaling function <p(t). Then the highpass coefficients produce the wavelets. 
Working with <p(t) and u;(f ), we really have three basic jobs: 

1. Compute the coefficients in fj(t) = J2k a Jk4>jk(t) and /(f) = T,jHk b jk w jk{t). 

2. Construct <f>(t) by actually solving the dilation equation. 

3. Connect the properties of 0(f) and w(t) to properties of the c's and <Fs. 
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Going between levels of a multiresolution is subband filtering with C T and D T : 
Theorems A function £a, +M * y+ u(0 * *«e space V j+i = Vj © W, has coefficients 
Oik and bjk in the new orthonormal basis {4>jk(t), Wj k (t)}: 

la; t = £c(£-2fc)a; + u ami &; t = £d(«-2*)a,+u- | (6.23) 

It f s 

/„ vector notaZ» thlsU a, = U 2) C «* »i - (1 2) D r ^ *y«»** * 





Proof. For j = 0, formula (6.23) is (6.21). The extension to every j comes from the dilation 
equation. Again n—l — 2k: 

2" 2 *(2>t -k) = 2" 2 £ *< 2/+l ' - 2* - n) = £>(*- 2*) 0, +w (t). (6.24) 

The wavelet equation has d in place of c. The inner products of these equations with /(t) give 
the recursions (6.23) for the coefficients a jk and b Jk . . 

Now go in the opposite direction. Change from the basis u>;*(0} back to the basis 

»/+W«). Since the bases are orthonormal, the inverse operation is given by the transpose. 

Theorem 6.3 a J+U comes from a jk and b jt by a synthesis fdter bank 

1 aj+u = f>(2* - 1) a ]k +d(2k -Qb jk J 
77te inverse pyramid is the fast inverse wavelet transform: 



b, 



(6.25) 



Lowpass Iteration and the Cascade Algorithm 

We begin the solution of the dilation equation. Our goal is to construct the scaling function <j> (t). 
The only inputs are the filter coefficients e(0) ..... c(JV). The first solution method we propose 

is the cascade algorithm. 

Start the cascade with *<°>(t) = box function on [0. 1]. Iterate the lowpass filter: 

^+i> (0 = £V2c(n)* (i) (2f -n) = J^2h(n)^(2t -n). (6.26) 

n n 

The algorithm woiks with functions in continuous time. Those functions are piecewise constant 
and the pieces become shorter (their length is 2"'). If^ (i >(0 converges suitably to a limit *(0. 
then this limit function solves the dilation equation. 
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Notice the two time scales, / and 2/ , which come from the continuous form of downsampling. 
In place of (I 2) <p(n) = <p(2n) y we have (4, 2) <f>(t) = 0(2/). The cascade algorithm is really 
iteration with the filter matrix Af = (| 2) 2JJ — as we will see in detail. It is an infinite iteration, 
and our final formula for <f> (/) will involve an infinite product 

It is easy to associate a continuous-time function *(/) with a discrete-time vector x(n). The 
function takes the value x(n) over the n th time interval. That is the interval n < t < n + 1. Thus 
the constant vector x = (. ..,1,1,1,...) produces the constant function x(t) == 1. The impulse 
x = (. . . , 0, 1, 0, . . .) produces the standard box junction. In general x(t) is piecewise constant: 
x(t) = x(n) on the interval n < t < n + 1. 

The iterations start from the box function 0 (O) (f). There are two steps in each iteration — 
filtering and reseating. Suppose the filter coefficients are h(0) = 2/3 and A(l) = 1 /3. Filtering 
the input gives 50 (O) (/) + 50 (O) (/ — *)• Then rescaling/ to 2t compresses the graph. Tomaintain 
a constant area we multiply the height by 2: 

* (,) W = f0 (O) (2O+|0 (O) (2/-l). 

Filtering and rescaling one box produces two half-width boxes of height f and §. That iter- 
ation step preserves the area (=1). Now filter and rescale 0 (1) (*)- The two half-boxes become 
four quarter-boxes, from 0 (2) (/) = f 0 (1) (2/) + f^ 0> (2/ - 1). The first quarter-box has height 
y . That height is multiplied by | at every iteration! 

We wish we could say that the iterations 0 (l) (*) are converging. Their limit <p(t) would sat- 
isfy the dilation equation </>(t) = f 0(2/) -f §0(2/ - 1). In some weak sense, this may be true. 
In a pointwise sense at / = 0, the functions 0 (,) (O) diverge because of (4/3)*. The coefficients 
2/3 and 1 /3 illustrate the iteration process, but not its convergence. 

We want to see that process also by algebra. It is clearest if we ignore the rescaling and just 
execute the filtering with coefficients h(k). The heights of the boxes would be f , A, and then |, 
| , | , | . In the z-domain, this corresponds to 

L H(z) = l + \z~ l and H(z 2 )H(z) - | + f z" 1 + -h jz" 3 . (6.27) 

|. The actual time intervals go from length 1 to \ to \. The actual graph heights are doubled at 
j each step, to preserve area. But the essential point is the product H(z 2 )H(z). After three steps, 
1 the iteration will produce ff 0) (z) = H(z 4 )H(z 2 )H(z). After i steps we have 

I 

^ ) W = flH(z 2t ). (6.28) 

>; This product is the z-domain equivalent of iterating the lowpass filter H(z). The values of 0 (O (/) 
fc— the heights of the graph after i iterations — are the coefficients of 2 l H (i) (z). That factor 2' 
J accounts for the height-doublings that preserve area, when the time intervals for 0 (,) (/) become 
S 2-' . 

- You may ask, why not choose the usual averaging filter as a first example? Let me show 
gg you why. The averaging coefficients are h(0) = h(l) = \. The first step of the iteration, with 
if coefficients 2A(0) = 2h(l) = 1, is 

%}' 

H 0 (1) (O = 0 (O) (2O + 0 (O) (2/ - 1). 

g From the box function 0 (O) (O this produces the same box: 0 (1) (O = 0 (O) (O- 

life 

K 
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The output equals the input The iteration process converges immediately. We have found 
the scaling function! In general 0 (0 is the limit of the sequence 0 (i) (O. when that limit exists 
as i -► oo. Here 0 (O) = 0 (1) and the box function is a "fixed point" of the iteration. When we 
filter and rescale 0(0 we get back 0(0, because the sum of two half-length boxes is the original 
box: 

(6.29) 



Box: 0(0 = 0(20 + 0(2/ - !)• | 
The z-domain equivalent is a product built from 

Please notice that we do not square this function. H (2) (z) is H(z 2 )H (z): 

H*ft> = (J + + = 1 + + + K 3 - 

After i iterations, HHz) will have 2< coefficients all equal to 7r*. After rescaling, this stiU 
corresponds to the box function. 

Now US e three filter coefficients h = (\,\,\). The box 0<®(O produces three half-boxes 

0<»>(O = l0< o >(2r) + 0 (O) (2r - D + \<t> (0) Qt - 2)- 
Then there are seven quarter-boxes in 0< 2) (O- Rescaling prevents the support interval from be- 
coming long. The limiting interval is 0 < t < 2. 



(6.30) 



•0 (,) (O 



0 {2) (O 





Figure 63: The cascade algorithm for \,\,\ converges to the hat function. 



A reasonable guess for the limiting function 0(0 is the hat function. This is piecewise linear, 
going up to 0 (1) = 1 and down to 0(2) = 0. We verify that the hat function is a fixed point of 
the iteration. Filtering and reseating leaves this scaling function 0(f) unchanged: 

0(0 = 10(20 + 0(2/ - 1) + ^0(2* - 2). (6-3D 

Notice how the coefficients \, \, \ aredoubled. The hat function is a combination of three nar- 
rower hats. For future reference, we note the different properties of these examples: 

1. H(z) = | + 5Z" 1 is not zero at z = -1, corresponding to o> = n. The iterations fail to 
converge. 

2 #( z ) _ I + I z -' is zero at z = -1. The iterations converge. The filter H (z)H(z -1 ) « 
" halfbandfno even powers except the constant term. The box function is orthogonal to its 
translates. 
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3. H(z) = \ + ±z' 1 + \z~ 2 is zero (twice) at z = -1. The iterations converge. The filter 
H (z)H(z~ l ) is not halfband. It contains the even powers z 2 and z~~ 2 . The hat function 
<t> (t) is not orthogonal to <f>(t - 1). 

We must quickly emphasize that a zero at z = -1 (which is co = n in the frequency domain) 
does not guarantee the convergence of 0 (O (f). But without that zero in the filter response, strong 
convergence has no chance. 

Similarly, a halfband filter does not guarantee that 4>{t) is orthogonal to its translates. But 
without that halfband property of H(z)H(z~ 1 ) i orthogonality has no chance. Section 7.2 will 
further indicate those connections; they are not quite two-way implications: 

Convergence of (t) to <p (t) needs J/ = 0atz = -l 
Orthogonality of <f>(t - k) needs H(z)H(z~ l ) to be halfband 

Orthogonal Functions from Orthogonal Filters 

When the filter bank is orthonormal in discrete time, we hope for orthogonal basis functions in 
continuous time. All wavelets w(2 j t) should be orthogonal to the scaling functions </>(t - *). 
Furthermore, the wavelets w(2 s t - k) should be mutually orthogonal and the scaling functions 
<}>(t-k) should be mutually orthogonal. Note that <f>(t) is not orthogonal to 0(20- 

Theorem 6.4 Assume that the cascade algorithm converges: # (l) (f) -> </>(t) uniformly in f. 
If the coefficients c{k) andd(k) come from an orthonormal filter bank, so they have double-shift 
orthogonality, then 

1. The scaling functions <p(t —n) are orthonormal to each other: 




Ht - n)<p(t -m)dt = S(m - n). 



2. The scaling junctions are orthogonal to the wavelets: 

f <p(t-m)w(t-n)dt = 0. 

3. The wavelets w jk (t) = 2^w(2h - k) at all scales are orthonormal: 

/oo 
u)jk(t)w JK (t)dt = 6(j - J)6(k - K). 
■oo 

Proof of 1: The box functions 0 (O) (t -k) are certainly orthonormal (because nonoverlapping). 
We will show that when <f>«\t - k) are orthogonal, the next iterates 0 (f+1 >(* - k) are also or- 
thornormal. Then the limits <p (t - k) are orthonormal. 

% The induction step from i to i + 1 assumes that the </> (/ >(f - it) are orthonormal, and sets 
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(6.32) 



= 2 J(£c(W (i) (2f - 2m - *)) £>(W (i >(2. - 2n - *)) dt 

= J(X>(W (i) (2.-2^^^ 

= £c(Jfc)c(fc -21) = 6(1) = 6(m-n). 

• t u 1act i; n e when we used the orthogonality of the row [c(0) • • • c(N)] 
Thecrucial step came m the last line, when we usea 8 . ^ statement 

to its double shifts. These are rows of L = 4 W ^^f nor ^ izeah alfband filter: 
LL r = / Equivalent^, it is in the statement that |£, c(fc)e | is a norauu 

the dilation equation would require = 0. 
Proof of 2: Repeat the integration steps above for <f> times w: 

J (f>(t-m)w(t-n)dt 

= /£>(*> V2 0(2. - 2m - k)) £>(*) V2 0(2f - 2« - *)) A 
= ... = J^c(k)d(k-2l)=0. 

„ The last steo uses the orthogonality of the rows of L = U 2)C to the 

^rltrixtbrraoftbisd™,^ 

„4e »(0 orthogonal to ♦(«)• To reach the end of part 2, we needed part 1 

between the ^>'s. , 

ProofofS: The orthogonality of wavelets ^0 at the same scale level (the same ,) is proved 

as in parts 1 and 2: 

/•"„(,_ „)„« -.)*—;- - m - m = m ' n> - 

Agato cooho^s ft. orftogt-* fo.lo»* *r£~"* ™* " 

2^ ia.orthogondfflteebarJ.mainta^^orthogo^ 
^tacUons will lead » the satrte ibted poht. *(.), or a. least to a trrttlaple *C0 
convergence holds. 
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In general, convergence can be "wea*" or "strong". For weak convergence, the functions 
<t> i0 (t) can oscillate faster and faster. You would not call this convergence. But the integral of 
converges to the integral of <p(t), on every fixed interval [0, T]. (Integration controls the 
oscillations.) In the convergence that we assumed, 0 (,) (O approaches <j>(t) at every point 

There is a better starting function <t> {0) (t) than the box. The constant value <p (0) (n) on each 
interval n < t < n -f 1 can be the correct <p(n). The values of <p are filled in at half-integers 
and quarter-integers by the iterations <f> il) (t) and 0 (2) (O- The graph of 0(0 appears, 2* points at 
a time. We stop when we have enough points for the printer to connect into a continuous graph. 

The next section explains how to start with the correct values of <p(n) at the integers. 



Problem Set 6.2 

L For the filter with h(Q) = h(l) = \ and any <f> (0) (t), describe and draw 0<'*>(f). 

2. If H(z) is a polynomial of degree iV, what is the degree of H(z 2 )H(z)l What is the degree of 

tf tf) (*)=nii 

Rescahng will replace z by z 1/2 . After/ steps, the degree is divided by 2\ Show that the degree 
of H {i) (z l/2 *) approaches N as i -» oo. 

3. With coefficients h(0), . . . , h(NX the support interval of 0 w (f) grows to [0, N]. What hap- 
pens if 0 (O) (0 is a box on [0, 2N]1 

4. The unit area of the box is preserved if and only if h(0) H +- h(N) = 1. Are negative 

coefficients allowed? 

5. Suppose the filter coefficients^*) are \, 0,0, \. Starting from the box function, take one step 
of the cascade algorithm and draw (0. Then take a second step and draw <f> a) (f). Describe 

| ; <t> {i) (t) — on what fraction of the interval [0, 3] does <£ (0 (t) = 1 ? 

6. Suppose the only filter coefficient is h(Q) = 1. Starting from the box function <p (0) (t), draw the 
graphs of 0 (I) (f) and 0®(f). In what sense does 4>&(t) converge to the delta function 5(f)? 
lb verify the dilation equation 5(0 = 25(2/), multiply by any smooth f(t) and compare the 

: g - integrals of both sides. 

f(5 7. Suppose 4> (0) (t) is a stretched box of unit area: ^(0 = 1/2 for 0 < / < 2. Draw the graphs 
of <f>V\t) and 4>™(t) when A(0) = h(l) = 1/2. On what interval is nonzero? What is 
the limit 4>{t)l 

■ ■ 8. Suppose 4>^(t) is the Haar wavelet with zero area: 

| <£ (0) (0 = 1 for 0 < t < 1/2 and tf (0) (0 = -1 for 1/2 < t < 1. 

. g Withfc(0) = = 1/2, draw the graphs of 0^(0 and>* 2 >(r). The sequence <j>^(t) con- 
verges "weakly" to what multiple c<f> (0? 

.*it<\t... 
>. - 

§•3 Computing the Scaling Function by Recursion 

; Jhe main point of this section can be stated in three sentences. Then you can follow through on 
^th| (letails, or look ahead for the matrices ro(0) and m(l): 

il ^ dilation equation gives #(0), <£(!), ... as the eigenvector of a matrix m(0). 
-Sp? ^ W at ' = half-integers comes from multiplying by a matrix m(l). 
: | -V Then <j>(t) at every dyadic t comes by recursion. Each step uses m(0) or m(l). 
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The scaling function is created recursively. This section gives the rule. 

The dilation equation is easiest with only two coefficients (N = 1). Then m(0) = 2ft(0) and 
m(l) = 2A(1) are scalars not matrices. The two-coefficient dUation equation is 

Itif) = m(O)0(2O + m(l)0(2t - 1). I t 6 - 3 *) 

The solution will be zero outside the interval 0 < t < 1. Inside that interval, set t = 0 to find 
0(O)andr = itofind0(|): 

0(O) = m(O)0(O) and 0 (\) = m(l)0(O). 

Now set t = i and J. then * and f , and onward through all the dyadic points t 
Directly from the equation you find 



= n/2 l . 



m 
♦ft) 
♦ft) 



m(Q)0ft) and 0(|) = 
m(O)0(i) and 0(f) = 



m(l)0($) 



= m(O)0(l) and 0(|) = m(l)0(f) . 



1 



Each new value comes from multiplying a previous value by m(0) or m(l). At each time f, the| 
right side of equation (6.33) has only one nonzero term. Thus 0(|) equals m(l)0( 5 ) which is|g 

"^iTtiletS step 0(f) equals m(0)m(l)in(W(0). The key is in the order of m(0) and m(ljj 
It is the sameorder as in the binary expansions f = 0.11 and f = 0.011. At any point* = n/2| 
the solution 0(0 has i factors: 



If f =0.01101 in base 2, then 



0(0 = m(O)m(l)m(l)m(O)ro(l)0(O) 

We have now solved the'two^oefficieat dilation equation at aU dyadic points 

Admittedly, the restriction to two coefficients looks severe. The pattern is correct and impop 
tant but two numbers m(0) and m(l) are not enough. The only normal case is m(0) = m(l) - v 
1, when we get the box function. Form(0) = § and m(l) = |, the first equation becomes,. 
^(0) = ±0(0). This produces a singularity of 0(f) at all dyadic points. /| 
We will not pursue that example here, because there is a more valuable application - which 
reduces N + 1 coefficients to two. This is the f amffiar step of reducing a high-order equatiojg 
to a low-order system. For differential equations that produces a matrix, asm the system « -g| 
Au For dilation equations the reduction will produce twomatrices m(0) andm(l). Theddatign| 
equation will become a two-coefficient matrix equation. The recursion will not change, exce|| 
it has vectors and matrices. '$£ 



Vector Form of the Dilation Equation 

The N + 1 coefficients in the dilation equation are V2c(*) = 2h(k): 



*=0 



Outside the interval 0 < t < N, we want and expect 0(0 S 0. Inside that interval, subjtig 
f = 0, t = 1, . . ., t = N - 1 to determine 0(0 at the integers. You will see again tne 



m 
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fact that even k goes with even 2t - k, and odd k goes with odd It - k. (Reason! The sum ofk 
andlt -k is even.) The right side of (6.34) leads to an N by N matrix m(0) which is displayed 
for N = 5: 



"0(0)" 




0(1) 




0(2) 


= 2 


0(3) 




0(4) 





MO) 

h(2) h(l) h(0) 
h(4) h(3) h(2) h(l) h(0) 
h(5) h(4) A(3) h(2) 
h(5) h(4) 



0(0) 
0(1) 
0(2) 
0(3) 
0(4) 



= m(0)$(0). 



(6.35) 



This is the dilation equation restricted to the integers. It is an eigenvalue problem for jb(0). That 
matrix has even and odd in separate columns. For a nontrivial solution, this matrix (including 
the factor 2) must have A. = 1 as an eigenvalue. Assume this is true. Then the values 4>(n) are in 
the eigenvector (which we call 4>(0)). That eigenvalue problem for m(0) sets the integer values 
0(n), and the recursion starts. 

Now look at the vector of half-integer values. Substitute t = \, §, §, \, \ into the dila- 
tion equation. This leads to a closely related matrix m(l). The first row comes from <j>(\) = 
2&(1)0(O) + 2ft(O)0(l). Notice that 2r is an odd integer, so the sum ofk andlt-kisnow odd. 
The matrix is m(l): 



" 0(1/2) " 




0(3/2) 




0(5/2) 


= 2 


0(7/2) 




0(9/2) 





h(l) &(0) 

h(3) h(2) A(l) A(0) 
h(S) A(4) A(3) h(2) *(D 
A(5) h(4) *(3) 
h(5) 



0(0) 
0(1) 
0(2) 
0(3) 
0(4) 



= m(l)<D(0). 



(6.36) 



As expected, the values at half-integers come from the values at integers. A vector <t»(i) comes 
from a vector *(0). In matrix notation (6.35) was an eigenproblem for <t>(0) and (6.36) is the 
step to $(5): 

4>(0) = m(0)«t>(0) and <&(I) = m(l)<I>(0). 

This is exactly like the two-coefficient case! Now m(0) and m(l) are N x N matrices. The 
beautiful fact is that the same pattern continues to quarter-integers and beyond. 

When t is a quarter-integer, the times 2t - k are half-integers. The values (f>(\), 0(f), ... 

come from <j>(\), (j>(\) The dilation equation connects those vectors by the matrix iw(0). 

Similarly the values 0(|), 0(2), . . . come from multiplying those half-integer values in the vec- 
tor O(i) by the matrix m(l): 

* (I) = OT (0) <&(!) and <*(!) = m(l)«t>(!) . 

Exactly as before, the binary expansion of t = n/2' reveals the order of the factors m(0) and 
m(l) — as they multiply the initial eigenvector 0(0) of values at the integers. We describe the 
recursion and then prove it is correct 



Theorem 6.5 The vector form of the dilation equation is 

4>(t) = m(0) <D(2f) + m(l) <D(2r - 1). 



(6.37) 
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J - o! 1 T • . . W - 1 are in the eigenvector of m(0) with X - 1: 
(Fixed point) *(0) = m(0)*(P). 

, j- ■ (+ w-l comes from i multiplications 

The vector <K0 of values at the dyadic points t, t+l t+N 

by m(0) and m(l). Here I = mflf < 1 with n odd: 

0(r + l) 



(6.38) 



Iff = |=0.011 then <D(0 = 



= m(0)m(l)m(l)<l>(0). (6.39} 



0(f + N-l) . 

sion, in which the 0-1 digit t x tells whether to use m(0) or m(l). 

V CCto rrec««ion: hh * . . .) = »(<l) •< •* 0 '4 -)• « M 

is nonzero for \ < t < 1. These two vec ^ fe substitute t into me dilatio 
are2A(2i-i): 

• ft(0) 

. ft(2) ft(D *(0) . (6.4 

. fc(4) ft(3) fc(2) *(1) *(0) 



Af = V2L = 2 



the dilation equation. 

Dilation Equation in Infinite Vector Form *«>(*) = M#oo(2t) 



This form is 



*(r-D 

0(0 
0(t + D 


= M 


0(2f - 1) 

0(20 
0(2f + 1) 









for - oo < t < oo(6.< 
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the vector form (6.37): 
*<0 = { 



m(0)*(2f) 
m(l)<J>(2r - 1) 



for 
for 



\<t<l j 



m(0)<I>(2r)+jn(l)<D(2r-l). 



The matrices m(0) and m(l) are N by N sections of the infinite matrix M For i 7-01 
- 1 the matrix entries are ' 

j m(0)y = My = 2*(2i - ;) and m(l) y = M ( , y _, = 2A(2i - j + 1). 



Proof. The verification is in three steps, first for M and then m(0) and then m{\). Ihopethe 

SSSSrt w y ,n , t0 r * e sum E a W(2/ - k) ■ ° r ^*(^> « m 

2)«<P(2f). We now follow each step: 

(Verify M) Row zero of ^(r) = Jlfd)^) is 0(r) = 2£ft(*W(2f - Jfe) 

Row one is <j>(f + 1) = 2£ h(k»Qt + 2 - k). The double shift works. 
(Verify m(0)) Restrict to 0 < r < I and keep only rows 0 1 N-l 

Sln l e *2?^L 1} x7c° ^ * (2/ + *° = ° We 0nly nC6d CoIunms 0 to /V - 1 of the infinite 
matrix M. This tV by N section is jb(0). 

(Verify m(l)) Restrict to i < t < 1 and keep only rows 0, 1 N-l 
Since <K2t - 2) = 0 and <P(2t + N - 1) = 0 we only need columns ' 
N — 2 of M. This /V by tV section is m(l). 

The change from m(0) to m(l) comes as f crosses I, because the nonzero entries in 
appear one component earlier. 

We now go back to the eigenvalue problem *(0) = m(0)*(0). Condition A\ leads to the 
eigenvalue A. = 1, and guarantees a solution. 

The Fixed Point Equation #(o) == m(o)$(o) 

The rowsof m(0) have a double shift. The columns have entirely even indices or entirely odd 
indices. The 5 by 5 matrix (N = 5) shows this pattern: 



m(0) = 2 



h(0) 

A(2) A(l) /»(0) 

A(4) A(3) A(2) A(l) A(0) 

MS) h(4) h(3) h(2) 

H5) h(4) 



^a^^T=T OD ^ C ° effiCientS Hin) iS C ° nditi0n A ' : The frequency response H(o>) 

A(0)-A(l)+A(2)-.. = 0. 
anT^st* A( ° ) + A(1) + * (2) + "' = thismeansthatev «y «*■» of mm andm(l) 



Theorem 6.6 Qwdftto, ,4, grantee, (teUl is an eigenvalue ofMandm(0) andm{\): 

*(0)-A(l) + A(2)-... = 0 ^ „ 

h(0) + h(l) + h(2) + ... = l Fi*2^A( B ) = 2j;i(fl) = i. 



odd n 



I 



llll 



- 1 -' -v 



mmm 



mm 
mm 



00 
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Any matrix with unit column sums has X = 1 as an eigenvalue. Theref ore ■. *&> = «(0)*(0) 

can ibe solved to give the scaling function at integer times $(0) = (0(O),0(1) K* ~ D) • 

Proof. Adding the two equations gives the even part Subtracting gives the odd part These 
aTL column sums of the matrix »(0), all equal to 1. A matrix wiffi unit column sum. has a 
left eigenvector of ones, e»(0) = e, because the multiplication just adds up each column. 

[ 1 1 ... 1 ]m(0) = [11 — . 1 ]• 

This means that«(0)-/is not invertible, and X = 1 is an eigenvalue. The left eigenvector is e. 
The right eigenvector is(0(O), 0(1), . . . , 0(M - W- 

(m(0) - /)*(0) = 0 which means *(0) = m(0)«I>(0). 

The fundamental fact is that a square matrix and its transpose have the same determinant and 
same rank and same eigenvalues. 

Tnecolumns ofm(l)andM also add to 1, producing the eigenvalue X = 1. SmaU point! We 
can^fety normalize the eigenvector 4>(0) by £ 0(n) = 1. This isthe "unit area" 
mat we impose on me mncnon 0<°> ^ 
f<(>(t)dt = 1 at the end. 

Corollary The sum^M + V h identically 1. 

Proof. Multiply the vector dilation equation *«) = m(0) *<?) + m(D *(2r - 1) on the 
left by e. Use the fact that e m(0) = e and e m(l) = e: 

e <p(t) = e 4>(2f) + e *(2* - !)■ 

This is a dilation equation for e <t> (0 and its solution is the box function! Thus e »(r) = 1 . The 
"periodized scaling function" £00 + k) is identically one. 

Example^. The coefficients ttflfc) = J. 1. J lead to the hat function. Tbe 2 by 2 eigenvalue 
problem for m(0) gives the correct values of 0(0) and 0(1): 

0(0) = 0 
0(1) = 1. 



[!■:][«]-[.«] - 



The sumofaUhatfunctions0(r-r«)is identically one. Noticematffiefirstrowofmee^ue 
equationisalways2fc(O)0(O) = 0(O).Then0(O)i S zoro,apartfrom 

^which occurs for the box function. This means that the scahng function 0 (0 - «o up to and 
including t = 0. The box function starts with a jump at t = 0, because fc(0) - 5. 

Example6.4. The Daubechies coefficients have 8A« = 1+V5.3+V5.3-V5.andl-V3- 
Dividing by 4 we have 2h(k), the numbers that enter m(0): 



1 + V3 0 0 
3-75 3 + V3 1 + V3 
0 1-V3 3-V3 





" 0(0) " 




" 0(0) " 




0(1) 




0(D 




0(2) . 




0(2) . 
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The eigenvector gives 0(0), 0(1), and 0(2): 

0(O) = O 0(l) = i(l+V3) 0(2) = i(l-V3). 

We now know the Daubechies scaling function 0(f) at the integers. The only nonzeros in 
the fixed-point eigenvector * are 0(1) and 0(2). From these values at t = 1 and t = 2 the 
recursion produces 0 (f ) at any dyadic point * 

Practical conclusion Every 0(/i/2') comes via m(0) and m(l). 

Theoretical conclusion Those dyadic values have a uniform bound if and only if all products 
ofm(0) andm(l) in all orders teve* uniform upper bound. When this holds, the dilation equa- 
tion has a bounded solution 0 (f ) for all f . 

I Wec^proposesufficientconditionssomataUproductso^^ 

| We can also propose necessary conditions. It is not known how to verify the necessary and suf- 
fficient condition (Section 7.3). 



tiff- 

Derivatives of the Dilation Equation 

|While working in the time domain, we might as well take the derivative of 0(f). The result is 
; : - highly interesting and not fully understood. Part of the problem is that the derivative 0'(f ) may 
not exist w 3 

:f f The plan is to differentiate each term in the dilation equation for 0(f): 

v m = 4£/H*)0' ( 2f-*). 

This is another dilation equation, with every coefficient doubled. The equation 4»(f ) = M \ 
has led to O'(f) = TMVQt). At , = 0 this yields the fixed-point equation *'(0) = W(0) 
J-The eigenvector 4>'(0) contains the derivatives 0'(n) at the integers t = n 
P|To solve <I>'(0) = W(0), we have a new requirement 77* number A. = i must also be an 
A ^ tWsa PP liestoth ^^-atrices W (0)and m( l ) . This new requu^m 
gKpe entnes is stated as Condition A 2 in the following theorem 

|pwrem6.7 The matrices M and m (0) andm{\) have eigenvalues 1 and \ if and only if the 
coefficients satisfy Condition A 2 which includes A,: 



<. Condition A 2 : jh(-l)'h(k)=0 and fV-l)**A(*) = 0. 

jf o o 

9f-"^ X = 1 * *?' C ° ntainin8 tk£ ValU£S *<°> *(» - D- The eigenvector 
• ~ 2 w<1> (OX containing the derivatives #'(0), . . - 1). 



»>is case //(a,) has a dmMe «, = „. This beautiful pattern extends onward to 

a^™5V^^r e eigCnValUeS ^ = M ^r'ifandonly if the filter 
Condition A p : £(-l)Vft(*) = 0 for m = 0, .. '., p - 1. 
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The eigenvector for X = (i) m contains values of the m* derivative of 0(0 at the integers. For- 
mally, <t> (m) (0) = 2 m M* (m) (0) conies from differentiating the dilation equation m times at the 
integers. We mention the frequency domain equivalent: 

Condition A p : The frequency response H{to) has a zero of order p at o) = it. 

We will see this again! And we also begin to uncover the crucial role of the left eigenvectors 
(which are row vectors). Those tell how to produce polynomials from combinations of the trans- 
lates 0(r - Jk). Under Condition A p , these low order polynomials 1 , t, .... t> 1 are in the low- 
pass space V 0 . They are the keys to approximation of a function /(f) by functions in V 0 . 

The letters A p indicate "approximation of order p. " The theorem above, with its extension 
from 2 to p, is absolutely basic to the algebra of downsampled filters. These eigenvalues and 
eigenvectors control everything in Chapter 7. 

You will see that the derivative 0'(O is often one-sided. Derivatives of 0(0 may not exist 
in the usual sense. This subject still contains some mysteries. 

Example 6.5. The hat function coefficients 2h(k) = 3,1,5 satisfy Condition A 2 : 

First sum rule: £ — 1 + 2=0 
Second sum rule: 0 ( 3) - 1(1) + 2 (*) = 0. 

Therefore m(0) will have eigenvalues 1 and \: 

, m f 2ft(0) 0 ] _ [ \ 0 ] 

m < 0)= [ 2fc(2> 2*0) _r|_ i 1 J" 

The eigenvector for X = 1 has components 0 and 1. They agree with 0(0) and 0(1), the hat 
function at the integers. The eigenvector for X = \ has components 1 and -1. They are 0^(0) 
and 4>' (1), the slopes 0'(f ) of the hat function in the two intervals. These are derivatives from 
the right at the points t = 0 and t = 1. The slopes on the left side of those points are different 
because the hat function has corners. 

The matrix m(l) must also have eigenvalues 1 and \: 



/is r 2A(i) 2A(0) 1 _ r 1 \ 1 

m(1) = [ 2*(2) J - L 0 5 J' 



The eigenvector for X = 1 has components 1 and 0. Those agree with the hat function at the 
shifted points t = 1 and t = 2. The eigenvector for X = \ has components 1 and -1. Those 
agree with the slopes of the hat function^rom the left at t = 1 and t = 2. Remember that m(0) 
is involved at the start of an interval and m(l) is involved at the end of an interval. 

Condition A3 is not satisfied for the hat function. There is no eigenvalue X = f The hat 
function has no second derivatives at the integers. 



Problem Set 6.3 



1. If the filter H(z) is halfband, show that the eigenvector in m(0)* = * is an impulse 6{n). 
What are the values of 0(t) at the integers? 
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' S^on? sadsfy equati0^s • does •* product = ^ « satisf y a 

• Show rthat the convolution^W^CO^satisfyadilation equation with coefficients from 
. Find a specific function /(r) that does not satisfy any dilation equation. 



Infinite Product Formula 

caling function 4>(t) comes from the dilation equation 

N 

Ht) = 2j2h(k)4>(2t-k). 
t=o 

f(t) I is the fixed point, or fixed function, when we iterate with H and rescale. In the time 
in, the matrix that filters and rescales is M = (4, 2)2* Now we intend to find the Fourier 
orat 0(a>) in the frequency domain. Just as the time-domain solution involved products of 
and m(l), the frequency-domain solution will involve an infinite product of H(co)'s 
ts quite remarkable that two-scale equations received so little attention for so long ' His- 

t ' \ u ? 0fteD S6en iD SamC eqU3ti0n - ^ to «W« Prominently 

u:tals which are self-sunilar. Now, multiple scales seem to be every where. We meet them 
' Jook through multirate filters- with two scales. Then the iteration leads to all scales 
the 2 s were removed, the dilation equation would be an ordinary difference equation 
□efficients are constant, so we look for pure exponential solutions When you make 
ibsutution, you are effectively taking the Fourier transform of the equation. The trans- 
urns difference equations and differential equations (and dilation equations) into algebraic 
ons. We do that now for the two-scale equation, and we watch how It leads to to 12 
» dilafon equation becomes ««) = (f ) ft f ). This leads recursively to an infinite 
:t for <t>{fo). This transform must be a sine function when h(0) = A(l) = 1/2 —because 
ae-domain solution Ht) is a box function. That sine function must be orthogonal to its 
anons by e , because the box function is orthogonal to its translates <j> (t - jfc) We 
d study orthogonality and also approximation, which is controlled by "zeros at n " These 
ties are now studied in the frequency domain. ' 

Condition O for orthogonality: 



\H(<o)\ 2 + \H(a> + x)\*~ m 1 in the frequency ^omau7 
I 2 T.KW{k-2t) = b\t) in the time domain I 



double-shift orthogonality of the lowpass filter coefficients. It connects to orthogonality 
scaling functions *(,-*). We use the word "connects" rather than "implies," because a 

condition is eventually needed to insure orthogonality in the limit of the iteration The 
am discrete time to continuous time seldom goes wrong, but it can 

is orthogonality will appear in Theorem 6.10 as a neat statement about the Fourier trans- 
ot</> and w: 

00 w 

£|£0 + 2™)| 2 = l and £ *(a> +. 2nn)w(a> 4- 2k n) = 0. 



202 



Chapter 6 Muttiresolution 



■ft 1 ... 



'.■~:y--'.'.-'- : . 



•••• • 

. - '-, .i, • 



1 



The other side of the theory is about approximation. This imposes a very different condition on 
thefc(*) and the polynomial H(o) = IX*^-' 4 ®. 



(1 4. e -i<B \p 
) Q(fi)). 



This factor (l + e^Y means that H{a>) has a zero of order p at a> = . We will prove thai 

this puts the polynomials 1, t f" 1 in the scaling subspace V„. They are combinations o1 

4>(f - n) and they are orthogonal to w(t - n). In the frequency domain, there is again a nea 
statement about the Fourier transforms of <f> and w: 

$ has a zero of order p at every a> = 2nn,n ^0 
S) has a zero of order p at zero frequency. 

The wavelet coefficients of a smooth function /(f) = £ b jk w Jk (t) decrease faster wha 
p is larger. The estimate is \b Jk \ = O (2^>). This is valuable for compression. This sectio: 
does the frequency-domain algebra, to solve the dilation equation and to explain Condition ( 
and Condition A p . 

Transform and Solution of the Dilation Equation 

To transform the dilation equation, multiply by «-**. Integrate with respect to t: 

r<t>(t)e-"»dt = 25>(*) f <t>{2t-k)e- iM dt. 

J-oo i=0 J - eo 

The left side is 4(a>). In the integral on the right, set « = 2t - k and t = (u + k)/2: 

f 0 - k)e-*"2dt = r me-^^du = e-W f (f ) . 

Instead of t and It, the transform involves o> and a>/2. The dilation equation becomes 

£(„) = (£fc(*)^ t/2 )?(f) = *(§)*(!) • 

This is the result of filtering and rescaling. Filtering multiplies £(<w) by H{a>). Rescali 
changes o> to a>/2. The scaling function (which is unique up to a constant multiple C—ti 
is still to be proved) comes out unchanged: 



(6.4 



(6.4 



&») = "(f) £(!)•! 



Now iterate this equation. It connects eo to a>/2 and therefore it connects a>/2 to <u/4: 

K») = *(!)["(!)*(!)]• 

After N iterations, this becomes 

?(a,) = H(f)H(f)-H(^)fe). 

In the limit as N oo, we have a formula for the solution ?(a>). Note that a>/2 N is approach 
zero, and 0(0) = f 4>(t) dt is the area under the graph of <f>0)- This equals one. We imp 
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the normalization 0(0) = 1 in the frequency domain, just as we required unit area in the time 
domain. Then the formal limit of the iteration leads to the famous infinite product for 0: 



(6.45) 



We note the minimum requirement for convergence of this or any infinite product: the factors 
H[a>/V) must approach 1 as ; -* oo. Thus we need H(0) = 1. By periodicity ff(2*) = 
H(4n) = 1. Then the equation *(a» = tf(f) ?(f) has a remarkable consequence. The 
values <K2jc), tf(4jr), 0(8jt), ...are all equal IfH(n) = 0, those values equal zero because 
<p(2n) = H(jc) <f>{n). 

This "zero at n" is a natural requirement on H(a>) in order that 0(<y) may decay and <b(t) 
may be a reasonable function. 

The infinite product converges for every a* and every H(a>). We have an explicit formula 
for 4>{a>). Whether any function <f>(t) has this Fourier transform is another matter! Convergence 
follows from a rough bound on H(a>) in terms of C = max | H'(co) |: 

\H(a>)\ = |1 +H(a>) - H(0)\ < 1 + C\a>\ < e CM . 
Then the product $(a>) has the same upper bound: 

M*>)\ = |f/(f )| |tf (f )| . . . < e c W 2 e c W* . . . = e CM . 



V This is a wild overestimate of $ (<»), as almost any example will show. 



tiff 



1 Box example. The coefficients are A (0) = ft(l) = I. Then H(a>) = Ul + e^) ^ 
? product of the first AT factors contains 2" terms. Looked at correctly, those terms are the fast 2" 
powers of e~ ta}/2 : " ~ ~ 



in-; 



H (m («.) = i (i+ ,-'»") (1 +.-**)... (i +,-**>») 



1 - e-' 4 ' 



2"-l 



C~4 • 

^KNbwIet^ 



2^(1-^/2^) 



, A 2 (geometric series) 
(sum of 2 N terms). 



(6.46) 



The denominator has = -/*+. . .) m w+ . . . with ^ = ^ 

: nie limit of 2«i0 is ico. Therefore the limit of the partial product is the infinite product 



Sftf '• 

Sill ^ J™* 5 * 011 is * e transform of the box function. The integral of from 0 to 1 agrees 

^ - " * * ^ me transform 



£fc>) = n (5 + |c- f<B/2 ') = (1 - e- to )/,a,. 



(6.47) 



.Sri;.' 



^p: actuaUy decreases to zero as co becomes large. 



g^^ompare the construction of 0(0 with In Section 6.2, we assumed that 0%) con . 

uniformly to *(r). Then we studied its properties. In this section, the convergence of 
ififej?&*te product is cheap (for each separate w). What we need is sufficient decay of ?(*>) as 
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\co\ oo. Our precise assumption will be continuity of the function A(co) in Theorem 6.10 be- 
low. Then we can safely study 0(a>) in the frequency domain. Note first that for real frequencies, 
the growth of 0 (a>) is at most polynomial: 

Theorem 6.8 |0(o>)l < e CN for complex co and |0(a))| < c(l + \a>\ M )for real co. 

Briefreason: H(co) is periodic. It has a maximum value 2 W . The equation 0(2a>) = H(co)<t>(co) 
says that 0 grows by at most 2 M when co is doubled. The bound \co\ M has this growth rate. A 
constant is included to make c(l + \co\ M ) correct for small |a>|. 

The example H(co) = \ + \e' m is bounded by 1 for real <o and by e M for complex co. Then 
M — o. The transform 4>(co) of the box function has those same bounds: 

^ I co I | co I [ 1 for real co 

|0(*>)| = |//(-)| \H(j)\ .-.<j e H/2 e M/4... = e M forcomplexo). 

Section 6.1 showed that the support interval for 0(0 is [0, N]. This can be proved in the 
frequency domain too. Our bounds on <j>(co) show two fundamental facts about 0(0: 

Theorem 6.9 Any dilation equation with h(0) + • • • + h(N) = 1 has a unique and compactly 
supported solution 0(f). This solution may be a distribution. 

Compact support comes from |0(<»)| < e CM . The Paley-Wiener Theorem implies that 0(0 
is supported on the interval [-C, C]. With more care [D, p.176] we could find again the exact 
support interval [0, N]. 

Uniqueness comes from our formula for the solution! The infinite product converges to 
0(a>), which is continuous because 0(0 has compact support: 

i oo ^ 

fii\co) = (fl H{co/2 j )) 0(a>/2 f ) approaches 0(*>) = (J\ H(co/2 j )j 0(0). 

' In the IIR case, suitable hypotheses will again give uniqueness (of course not compact support). 
At the other extreme, note how the lazy filter with A(0) = 1 leads to 0(0 = delta function. The 
dilation equation 0(0 = 20(20 is solved by 0(0 = 3(0: 

In frequency: H(co) = 1 so 0(a)) = 1. 

Cascade algorithm: 0 (O (O = box function on [0, 2~*] with height 2\ 
Verify directly: 8(t) = 25(20 from / f(t)S(t)dt = /(0) = / f(t)S(2t)2dt. 

All these methods show that ft(0) = 1 produces the best-known distribution 0(0 = S(0- 
Orthogonality in the Frequency Domain 

The product formula for 0(a>) applies with or without Condition O. When that condition holds, 
we expect orthogonality of the translates <f>(J - k). To establish this orthogonality in the fre- 
quency domain, we need to know that the equivalent statement is A(co) s 1. The function A (co) 
enters naturally into this discussion. It is the transform JX*)^' 0 * of the vector of inner prod- 
ucts of 0(0 with 0(f - k): 
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Theorem 6.10 The inner products a{k) are the Fourier coefficients of the In-periodic Junc- 
tion A(o>): 

/<» oo 
<f>(t)(t>(t-k)dt transforms to A(a>) = V) \$(a) + 2nn)\ 2 . 

The translates 4>{t -k)are orthonormal if and only ifA(o>) = 1. 

Proof. An inner product in the time domain equals an inner product in the frequency domain, 
by ParsevaTs identity. The inner product in the time domain is between <p(t) and <p(t - *). The 
transforms of these functions are $(o>) and e- ieok $(co). Each inner product integrates one func- 
tion times the complex conjugate of the other: 

[ 4(*)M-k)dt = [™$(<o)4>We io *d(o 

J-oo 2Jt J-oo 

= ^j[ te f)|?(« + 2]r»)|V-*rf». (6.48) 

Hie last integral split (-oo, oo) into an infinite number of 2* -pieces, using the periodicity of 
e iak . This integral defines the Fourier coefficient of A(co). Thus A(co) = J^a(k)e ia,k . 

For an FIR filter, <f>(t) = 0 outside the interval [0, W]. The inner products are a(fc) = 0 
for |A:| > N, because 0(0 and 0(f - *) have no overlap. The function A(o>) = JX*)*** 
w a trigonometric polynomial of degree N, which is not obvious from £ |#(a) + Inn)] 2 . In 
Section 7.3 we will compute a(k) direcdy from the coefficients h(n). TTiis is always a main 
point of the theory, to return every calculation to those numbers h(ri). 

When the translates are orthonormal, all inner products a(k) are zero except for a(0) = 1. 
The function with those coefficients is the constant function A(co) == 1: 




We now apply Condition O in the frequency domain to deduce this orthogonality of <f> (t -k). 
We are repeating in the frequency domain the result of Section 6.2 in the time domain. I believe 
this is worthwhile! The arguments in the two domains look quite different Recall the condition 
on the frequency response H(co) to produce an orthonormal filter bank: 

Condition O: \H(a>)\ 2 + \H(co -h n)\ 2 = 1. 

This function H(co) leads to which leads to A(a>). Somehow, Condition O must imply that 
A(co) = 1. The steps are typical of computations in the frequency domain. 

Theorem 6.11 If A(o>) is continuous, A(a>) = 1 is equivalent to Condition O. 
Proof. We use a very important two-scale identity, proved below: 

A(2a>) = \H(a>)\ 2 A(co)+ \H(co + n)\ 2 A(co + n). (6.49) 
If A(a>) m 1, this immediately gives that \H (a>)\ 2 + \H(a> + n)\ 2 si. 
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For the converse, suppose that Condition O holds. The identity says that A(2a>) is a weighted 
average of A(co) and A(co + n). At the point (Oq where A(2a>) reaches its maximum, A(o>o) 
and A(coq + n) must also reach that maximum. Now repeat the argument at a)o/2, to show that 
A(aW2) shares this same maximum with A(cu 0 ). Continuing, the maximum of A(o>) is achieved 
at a>o/4 and coo/i and eventually (by continuity) at co = 0. 

By a similar argument, which is due to Tchamitchian, the minimum of A(a>) is also attained 
at Q) = 0. Therefore A(o)) is constant We verify below that the constant is one: A(o>) s 1. It 
only remains to prove (6.49). ^ 

This valuable identity for A(2a>) has a nice proof. It uses the dilation equation <p(2a>) = 
H(<0)4>{(0). At the points 2a> + Inn, this splits into separate cases for even n and odd n: 

$(2<o + 2nn) = H(o) + nn)^(co + nn) 

| ff(a>)£(G> + 2*7T) ifn = 2* 

~ { ff(a> + *)£(<» + (2* + l)?r) ifn = 2* + l. 

Now square both sides. Sum from -oo to oo on n and therefore on k. The sum of squares is our 
function A(2<a>) in the desired identity: 

A(2o>) = \H(a>)\ 2 J $(a> + 2nk) | 2 + \H(a> + n)\ 2 J \$(co + n + 2nk)\ 2 

= \H(a>)\ 2 A(co) + \H(a> + n)\ 2 A(o> + n). (6.50) 

The final step is to confirm that A(0) = 1. This comes from our other condition on the 
lowpass filter, not yet used in the frequency domain. Condition A\ is H (n) = 0. In the time 
domain, this first sum rule guaranteed an eigenvalue k = 1 for the matrices M and m(0) and 
m(l). The fixed-point equation 0 (1) (n) = 0 (O) (*) at the integers could be solved. Condition A x 
is equally essential in the frequency domain. Here we use it to pin down the value A(0) = 1. 

Theorem 6.12 // H (n) = 0 then 0(2jrn) = Ofor all n ^ 0. Therefore 

A(0) = £ $(2nn)f = |£(0)| 2 = 1. (6.51) 

-oo 

Proof. The infinite product for 0(2tt) = H(n)H(n/2) • • • starts with the factor H{n). Imme- 
diately this product is zero. For any higher value n > 1, write n = 2 j m with odd m. Then the 
( j _|_ i)« factor in the infinite product is zero when (o = 2nn: 

K2„n) = H{nn)H^) • • • HQ) ■ • • = 0 

because H(nnfV) = H (jrm). By periodicity this is H(n) = 0. The only nonzero term is 
|£(0)| 2 . But£(0) = /J(0)//(0)ff(0) • • • which is 1. 

Orthogonalization of the Basis 

The condition for an orthonormal basis is A(a>) = 1. When this is not satisfied, there is an easy 
way to make it satisfied. In other words: when the translates <p(t-n) are not orthonormal, there 
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is an easy way to make them orthonormal. Divide <t>(a>) by the given A((o) (or rather, its square 
root) to get the new orthogonalized function ^orthC^) 1 



^orthfa) = 



This immediately gives orthogonality of the new basis {<£ O rth0 — n )}- 



^orthM = £ 



<p(co + 2nn) 



M«>) 



si. 



That succeeds if A(a)) is never zero. This is the condition for & Riesz basis. 

Theorem 6.13 The upper and lower bounds on A(co) are the Riesz constants B and A for the 
basis {<f>(t - k)} of V 0 . Thus A(a>) > A > 0 gives a stable basis, and dividing #(o>) by ^/A^o) 
gives an orthonormal basis. 

When 4>(f) comes from a dilation equation — this is our normal situation — Condition E in 
Chapter 7 gives an equivalent test for a Riesz basis (in terms of eigenvalues). If this test is passed, 
the wavelets w jk (t) are a Riesz basis for L 2 (R). 

Proof. To test the linear independence of the functions <f>(f - it), form the matrix A from their 
inner products. The entries are Ay =' (<f>(t - i), <f>(t - j)) m That number is a(j - /), and A 
is a Toeplitz matrix! It is the matrix TT* in Section 2.5. In the frequency domain it becomes 
multiplication by A(co). The upper and lower bounds on A((o) determine whether {<f>(t - k)} is 
a Riesz basis. 

Orthogonalization is always a basic step in linear algebra There it is done by the Gram- 
Schmidt algorithm. We start with independent vectors and produce orthonormal vectors (or 
functions). This algorithm is not successful here, because it is not time-invariant. The orthog- 
onalized functions will certainly not be translates — when the Gram-Schmidt algorithm works 
on functions in a definite order like <p(t), <f> (t — 1), 0(f + 1), . . . . To keep a shift-invariant basis, 
y/e needed to orthogonalize all these translates at once. The division by ^/A(co) did it 
H In matrix language, M^Af^ = /. In the improved factorization by Fourier methods, all 
rows of Af ortn come from the zeroth row by double shifts. In other words, M ortn comes from a 
filter. 

\ One problem with dividing by JA(a>). This destroys the finite response of the original filter 
fl:^The ne w filte r JI ort h is HR, not FIR. The new scaling function 0 O rth(O that corresponds 
S> 0(a>)/VA(a>) does not have compact support Vetterli and Herley noticed that this is not 
as it see ms. Since <f>(t) is zero outside the interval [0, N] 9 the inner products a(k) = 
t)$(t - k) dt are zero for \k\ > N. The function A(o)) with these Fourier coefficients is a 
ireal non-negative trigonometric polynomial of degree N. Its square root G(co) = Yls(b) e ~ ika> f 
spectral factorization, is also a polynomial of degree N. Tht frequency response of the new 
pnalized filter is a ratio of polynomials 




is; 

m 



#orth<» = 



G(a>Y 



|||| input-output equation j>(*) = £ft 0 rth - *) is an implicit difference equation, from 
an autoregressive moving average filter: 



,0 



■mm 



If 



m 



m 
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The new filter is HR but it only involves 2N +2 parameters g(k) and h(k). Therefore it can be 
physically realized, and now the basis has been orthogonalized. 

Problem Set 6.4 

1. Use the identity sin 29 = 2 sin 6 cos 6 to show that 

/ co\ / Q)\ ( o>\ 1 sina> sin f sin 
(co Sl )(cos-)-(co S ^j = ^ — 

Cancel sines and let N oo to find a great infinite product 

n c ° s (5) = ^ sinG> - 

2. The Haar filter has H(a>) = |(1 + e~ io> ) = e'^ 2 cos f . Use f in Problem 1 to give a new 
proof for the infinite product (6.47) of H(a>/2 J ): 

(«--|)-=(«-'-»!^?)=i(—-)- 

3. Suppose H(co) = }(1 + <r ,4> ) 2 . Find £(a>) and 0(f). 

4. If H(o>) has p zeros at o> = n> show that #(a>) has p zeros at a> = 2nn for each n ^ 0. 



6.5 Biorthogonal Wavelets 

This chapter has concentrated on orthogonal wavelets, coming from an orthogonal filter bank. 
The synthesis filters are transposes of the analysis filters. One multiresolution is all we need. 
The synthesis wavelets are the same, in this self-orthogonal case, as the analysis wavelets. But 
from biorthogonal filters we must expect biorthogonal wavelets. 

We now meet a new scaling function £(*). Its translates £(r - k) will span a new lowpass 
space Vb— different from V 0 - There is also a wavelet 2(0- The translatesJS(r - k) span a 
complementary highpass space W 0 . The sum of those spaces will be V x = V 0 + W 0 , the next 
space in the second multiresolution. To a large extent, the theory is achieved by inserting a tilde 
where appropriate. We want to indicate why this second scale of spaces Vj is needed. 

Biorthogonality comes automatically with inverse matrices. The rows of a 2 x 2 matrix and 
the columns of its inverse are biorthogonal: 

f row 1 1 f column column 1 _ T 1 0 1 
L row 2 \ [ 1 2 J-[ 0 1 J" 

Notice something pleasant The product in the other order is still /. The right-inverse is also the 
left-inverse. This order involves columns times rows, which are full matrices: 

j" column J [ row t ] + |- column j [ row 2 ]==L 

Those two column-row products are projections, and they add to /. These simple facts about 
2x2 matrices have important parallels for biorthogonal filters and biorthogonal wavelets. 
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Filter banks display those parallels immediately. The analysis bank has filters (| 2)H 0 and 
U 2)J?i. Those are rows 1 and 2. The synthesis bank has expanders before filters, F 0 (t 2) 
and Fi (f 2). Those are columns 1 and 2. Li the orthogonal case F 0 is #J and Fi is H T . In 
the biorthogonal case we don't have transposes but we still have inverses. To understand the 
pattern for wavelets, we absolutely must return to multiresolution. One scale of spaces V Q c 
V\ C • • • C Vj is too limited. We need two hierarchies of spaces, Vj in synthesis and V.- in 
analysis. 

Hide Notation Does the tilde go on the analysis functions or the synthesis functions? Both 
conventions are equally possible. We hope to agree with other authors! More and more, the 
tilde is going on the analysis functions. Then f(t) is expanded in synthesis functions, which 
have no tilde. But the coefficients come from the analysis functions and have tildes: 

fo(t) = £ 2o* <l> (t - k) is in V 0 , with = ff(t) $(t - Jfc) dt (6.52) 
f(f) = £ £ bji wjk is in L 2 , with b jk = //(/) w jk (t) dt (6.53) 

What does this mean for the filter banks that process the coefficients? Those filters use the letter 
H in analysis and F in synthesis. We will stay with H and F (rather than C and D) when dis- 
cussing biorthogonal filters. An important result in this section is the Fast Wavelet Transform in 
equation (6.70), and its inverse (the biorthogonal IFWD in Theorem 6.16. 



I Biorthogonal Multiresolution 

I This chapter began with orthogonal bases {0(f - *)} for V„ and {w(t -Jfc)} for W 0 . The equation 
| Vb 0 W 0 = Vi started a multiresolution. Wj was the orthogonal (!) complement of Vj inside 
Vj+i . All is well if <j>(t) and w(t) come from an orthogonal bank of FIR filters. Their translates 

an* all rvrtVirvrrrvTiol TtiAn J: 1 3 * 



are all orthogonal. They span perpendicular spaces and we have an orthogonal multiresolution. 
h AU is not so well if 4>(t) and u>(r) fail to have compact support. The filters fail to be FIr! 
JDften this means mat we have asked for too much! Instead of orthogonal bases, we should be 
g content with stable bases. An outstanding example is the space of piecewise linear functions. 
1 The stable basis consists of the hat function <f>(t) and its translates. That basis is not orthogonal. 



■ - - . ■ ..~ 6 ~.*u... 

I When the basis is not orthogonal, there is no reason to insist that W 0 must be orthogonal to 
V °- u we *>• m e multiresolution is called semi-orthogonal in Section 7.4, and we have "pre- 
Igwavelets". But the important property is a stable basis {w(t - Jfc)}. The highpass coefficients 
| will construct w(t). 

| H Remember the pattern for perfect reconstruction. Coefficients are chosen so that F 0 (z) H 0 (z) 
|,is halfband. When <f>(t) is the hat function from F 0 (z) = (^) 2 , the other factor H 0 (z) needs 
£p. coefficients. This means N = 2 but N = 4. The wavelet has 3-interval support. Then 



® md w( - f ~ s P m v o ^ Wo, without orthogonality 
■ | The new analysis multiresolution is the point of this section. The coefficients from H 0 (z) 
go into a different dilation equation, whose solution is the analyzing function £(i): 



N 

Analysis Dilation Equation: 0(f) = 2fto(*) ?(2f - jfc). 



(6.54) 



^^|)efficients Ao(A:) add to 1 as before. The new multiresolution obeys the same conditions as 
Wm^' m add a mt - V ° is SP 8011 ^ by [$(t - k)}. The space V, is spanned by {£(2'f - Jfc)} 



ill 



in 



fill 
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They are clearly shift-invariant. The dilation equation (6.54) says that V 0 C V,. ITienalso V) C 
V,+i . The highpass coefficients produce the wavelet: 



Analysis Wavelet Equation: w(t) = £ 2fci (fc) £(2f - *) 



(6.55) 



This wavelet is supported on [0, fl. where 21 = N + N. That sum N + is the degree of 
the product filters F 0 (z)H 0 (z) and *<*>«<*). Ttao are symmetric halfband filters and m the 
hatLction example thedegree is N + N = 2+4. Then* = Sisodd. The four functxons^t), 
u>(0, £(0> S(0 are graphed in Figure 6.6. 



Analysis Scaling Function 



Synthesis Scaling Function 





Figure 6.6: Biorthogonal scaling functions and wavelets from a 5/3 PR filter bank. 



Biorthogonality in Continuous Time 

Ourconstractionof^O^ 
coefficients*^ are notdouble-sWfto 

onal to the lowpass synthesis coefficients / 0 (fc): 

2 X>o(Wo(* + 2n) = 6(n). < 6 * 

This means that F 0 (z) flbfc) is halfband. Similarly the highpass filters give F, (z) ff, (z) = ^ 

band: si 

2][>(W,(* + 2n) = fi(n). 
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The other key relation is biorthogonality of highpass to lowpass: 

5>o(*yi(* + 2n) = 0 ^ £A,(jfcy 0 (jfc + 2n) = 0. (6.58) 
The reader knows that all these equations restate perfect reconstruction: 

F 0 (z) = ffi(-z), F,(z) = -flb(-z), F 0 (z)H 0 (z) is a halfband filter. 

Our question is, how does biorthogonality appear in continuous time? The functions #(r) 
and <j>(t) come from iterating the lowpass filters F 0 andfT 0 (with reseating!). Start the cascade of 
iterations from 0<°>(r) = £<°>(r) = box function on [0, 1]. Their translates have biorthogonality. 
The box 0<°>(r-*) has no overlap with £«»(*-*) when k£l. This biorthogonality is preserved 
at every iteration step, when we use equation (6.56). This is exactly parallel to the earlier proof 
(6.25) that orthogonality is preserved at each iteration. When ^<'>(f) and £<'>(r) converge in L 2 
to the scaling functions 0(f) and £(f), those limit functions inherit the same biorthogonality: 



f 

J—c 



4>(t - k)$(t - i)dt = 6(k - I). 



(6.59) 



Withi -> oom me cascade algorithm mese limits 0^^^ 
dilation equations. Now bring in the wavelet equations: 



f_J<f)${t)dt = jH (]T> o (*)0(2r -*)) (£2/i(W(2/ - €)) dt. 



(6.60) 



That right side is zero because of (6.58) and (6.59). And biorthogonality extends to the translates 
for the same reason: 



/oo 
•CO ' 



4>(t-k)w(t-e)dt = 0 for all* and £. 



: - Fmally the wavelets u» and 55 are biorthogonal from the wavelet equations and (6.57): 



r 

J— oo 



w(r-*)5(r-€)rff = $(*-£). 



(6.61) 



(6.62) 



All this is routine, provided the cascade algorithms for <f>(t) and$(t) both converge in L 2 . Wave- 
i*S g meor y 8> ves requirements in Section 7.2, as tests on eigenvalues of two matrices T and 
^^Suppose those tests are passed (not at all automatic!). The basis functions are biorthogonal. 
"j 1 ** does this say about the subspaces they span? 



Theorem 6.14 



Suppose the filter coefficients satisfy (6.56H658) and also Condition E (for 



^ik? 0nVersence °f the cascade algorithm). Then the synthesis functions <P(t - k), w(t - k) 
jgkfff^onalto the analysis functions $(t -£), w(t - t)asin (6.59H6.62). Each scaling space 



are 



t jj^hogonal to the dual wavelet space: 

Vj ± Wj and Wj ± V). 



(6.63) 



Jl|fe d ^ ^ P^endicular because their bases^a-Jfc) and S5(f _ k) are perpendicular. When 
^^^placed by Vt, the zero inner products are still zero. (Change variables back to T = Vt.) 
if Ifel-^ ^ SCaIe 7 we have P er P e ndicular subspaces. 
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W, 




Figure 6.7: V, is perpendicular to , while Vj is perpendicular to W). 
The two multiresolutions are intertwining. As always we have 

Vj + W J = Vj+ l and Vj + Wj-Vj+i. ' 



(6.64 



These are direct sums but generally not orthogonal sums. The subspaces V) and Wj have zer 
intersection, but they are not perpendicular. Instead Vj is perpendicular to Wj. All the subspaa 
W)-i, W)_2. - . . are then perpendicular to Similarly all the subspaces Wj-u W)-2t . . . ai 
perpendicular to Wj. Therefore we have biorthogonal bases (dual bases): 

Corollary The wavelets w jk (t) = V' 2 w{2h-k)andw jk (t) = 2>' 2 w{2h-k)arebiorihoi 
onal bases for L 2 : 

r w jk (f)w JK (t) dt = 6(j - J)S(k - K). (6.6: 

J— CO 



Representing f(t) in a Wavelet Series 

If we have a wavelet basis, we have a wavelet series. Any square-integrable (finite energy) fun 
tion /(f) can be expanded in wavelets: 



(6.6 



The synthesis wavelets are used to synthesize the function (of course). But the coefficients I 
come from inner products with the analysis wavelets. This is why b jk hasa tilde. Multiply (6.t 
by the analysis wavelet wj K (t) and integrate over t. Biorthogonality yields 



bjK = / f(f)w JK {t)dt. 
J -co 



(6.C 



Equations (6.66) and (6.67) are the biorthogonal wavelet transform and its inverse. The tra 
form takes function to coefficients, the inverse transform synthesizes the function. We show I 
low how the coefficients can be computed recursively (or pyramidally). This is the fast wave 
transform. 
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""^^'s^^ecnf l/«l 2 and EE|^| 2 is m>r^. ^required 
orthogonality - and not ^orthogonality. When we square and integrate the series for f(»Zn- 
zero inner products come fiom the products w jk {t)w JK (t). We do have an inequality 

A jy^ 2d ^EE\h\ 2 <Bj_y^^ (6 . 68) 

With B >A > Othissays we have a stable basis or /Km W This is true for wavelets subject 
to the same Condition E [Cohen-Daubechies]. waveiets, subject 



Fast Biorthogonal Wavelet Transform 



■W; 



The reader knows that in practice the wavelet expansion cannot go all the way back to / = -co 
WedonotusearbitrarUylowfrequencie^ 

sion starts wim V 0 and V 0 at a scale normalized to A/ = 1, and goes to V> and V>, wh"e Ae 
/ rt ED °^ h I hi ^ uenc y - included to reproduce an accumte sig^T 
Thus we work primarily with the subspaces Vj = V 0 + W 0 + ... + Wj , There are two 
nnportantbasesfor Vj . Oneis^(r) = 2^t - k) for -co < * < Jo T^e^ 
fimctmns are or /eve/ /. The other basis consists of M ) from V 0 and Wjl( t) for 0 < / <7 
Smce . kfe . recursive, we are interested first of all in / = 1. Then the two ways to'exoand a 
S1 gnalm V, are the scaling basis (fine scale) or scaling functions plus wavelets (coarse La^e): 

£ ^*<* " *> = X> ~ *) + £ *<* - *). (6.69) 
The key is me pyramid algorithm. This connects a lt to a* and We are using * and u, be- 



2^ = / f(t)$(t-k)dt= f /(0£Ao(<-2*)£ u (,)A. 

™ 1116 Wa ? l6t 6(1,13110,1 ^ *'« - a > instead «f the dilation equation with 
^ no(i - 2£). The Dvram H has fi1ws nA ^ ^ 1 1 



B ■ 

mi 

;; Ao(£ - 2*). The pyramid has filtering and downsampling- 

t§lW t^ t ^ ekt ^^ m ^k = Eho(e-2k)a lt and ^ = F Al (£ ^mz', 7% (67m 
; r :/^«^^~,«^JL±"' 1 2 * ,fl "' fl (6 - 70 > 



L-;- . - — - - - ; 

.Thisincludes a time-reversal! The same filters -fl£ and £Tf operate atleVel} 



- TV, of u ir. • « 0 upcraie ai level y. 

* ' pose^e^ r ^ 0rth ° g0nal - G ° ing backwards ' me »vme will not be the trans- 
P 086 - ™ e synthesis filters F 0 and F, must do their part 



,:\ - ---- * ' 

I V ^erseFast wavelet Transform: The coefficients a u in the basis M t) can 

computed from a^andbotby time-reversed synthesis filters: W 

git-Et/otf- tt)^ + E t /,(€-2*)^. I (6.71) 



jSlffp Perfect reconstruction operates when (6.70) is substituted into (6.71): 

ZU = " ^ - ^ + I> - I>« - 2.*. (6.72) 
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H^b, H^b/F, b.-^F, 

Figure 6.8: The pyramid algorithm from H 2 back to ai. 

The double-shift biorthogonality in (636) and (6 .57) makes this correct. The beauty of this Mai- 
lat algorithm (Figure 6.8 goes up and back down) is the way it connects continuous-time mul- 
tiresolution to discrete-time filters. 

Notice that the whole pyramid operates equally well if F 0 and F x are exchanged with H 0 and 
H\. The dual expansion of f(t) is: 

/(0 = £X>t £,*(') an d b jk = J f(t)w jk (t)dt. (6.73) 

All products have a tilde multiplying a non-tilde! This starts with the inverse relation of synthesis 
to analysis. Tildes can be exchanged with non-tildes throughout (if we want to do it). We will 
select <f> and w to be effective in synthesis. 

We emphasize one final point, often ignored. The coefficients in the expansion of f(t) are 
really differentfromsompfcs of /(*). They are inner products, not point values. This distinction 
is made in Section 7.1. 

Filter Construction by Lifting 

Herley and S weldens have proposed (independently) a systematic way to construct biorthogonal 
filter banks. Only one lowpass filter changes at each step. Starting from short filters, he quickly 
builds longer ones. In all cases the highpass filters remove aliasing in the standard way: ffi (z) = 
Fo(-z) and Fi(z) = -Ho(-z)- Then equation (4.9) on H 0 and F 0 is the remaining condition 
for perfect reconstruction. We drop the subscript zero on these lowpass filters, and recall the 
condition (4.9) that removes distortion: 

F(z)tf(z)-F(-z)ff(-z) = 2z-' (oddi). (6-74) 

Suppose this is satisfied by F and H; the filter bank is PR. Keeping F fixed, what are the other 
possible choices for H ? The answer is simple and important: 

Theorem 6.16 (Lifting) For fixed F(z), the solutions H*(z) to (6.74) are 

S H*(z) = H(z) + F(-z)S(z 2 ) for any S(z). j (6-75) 

Proof. Substitute ff # (z) to show that equation (6.74) is still satisfied. The new terms are 
F( z )F{-z)S(.z 2 ) ~ F(-z)F(z)S(z 2 ) = 0. This is in [HerVet] and [Swel]. 

Note that with F fixed, the equation is linear in H. We are starting from a particular solution 
(right side = 2z~'). To this particular H we are adding solutions F(-z)S(z 2 ) to the homoge- 
neous equation (right side = zero). Thus the even S(z 2 ) displays the degrees of freedom that 
remain when the PR condition is satisfied (Problem 5). That freedom is used in the Daubechies 
construction to achieve zeros at z = — 1. 
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We still want and need zeros at -1, which is <o = * . We also need stable bases. ( Section 
7.2 will state the stability requirement as Condition E. There is not yet a simple way to decide 
which S(z ) are permitted by this condition. In practice, we construct a potentially useful H*(z) 
and test it for Condition E.) This section will build in other important properties - linear phase 
interpolating scaling functions, binary (dyadic) filter coefficients. Those come at the expense of 
higher-order zeros at n. 

Dual lifting is also useful. In this case we fix the analysis filter H(z). The PR condition 
(6.74) becomes linear in F(z). The lifted solutions are 

F # (z) = F(z) + H(-z)T(z 2 ) for any 



T(z). ■ (6.76) 
Starting from the Haar filter or even from the 'Lazy filter", we alternate lifting and dual lifting 
to construct high-order biorthogonal filter banks with good properties. All these filters would be 
attainable directly from the (second-degree) PR equations. Lifting is a way to solve them as a 
sequence oflinear equations, with For H fixed at each step. Then we control more closely the 
final result 3 

Sweldens also emphasizes that Ufnngyiel(isa/art*ri>^^ 
! and its inverse. The Mallatfilter tree, which is subband filtering, is climbed in smaller steps This 
I is related to a lattice factorization. 

^ Example 6.6. The Lazy filter has H(z) = 1 and F(z) = z">. There is no true filtering, only 
g subsamplingfromU 2). Section 4.3 displayed block diagrams of this filter bank. Its polyphase 
^..matrix is H D = /. 



Suppose we keep H(z) = 1 and apply dual lifting to the synthesis filter: 

F*(z) = z- l + T(z 2 ). 



|.F*(z) can be any halfband filter, with one odd power z"' _ „ . ^ , u „ , 

contain powers of z as well as «-». This is needed in order to create symmetric filters. 



(6.77) 



z~ l . We are allowing S and T to 

•=:•. c , » create symmetric fihm. 

^Earher we centered the product H{z)F®, multiplying by z>. In this case centering gives 

i life 7« + f (2> - Th6n Da C ° meS fr ° m - C"« + 9 + 9 *~ l - * 2 )/16. Every maxflat 
gj|jgjfend filter D 2p can be lifted and centered from the Lazy filter. This section will create sym- 
metric biorthogonal filters, by lifting H = 1 while F = £> 2p . 

~ * 0te * e ! c ^S^^^ In analysis we have the delta function 

|||= 5(f) from H(z) = 1. This solves the dilation equation 5(r) = 25(2/). In synthesis F* 

yields an interpolating scaling function, with *(n) = 6(n). This is certainly biorthogonal to 
/ the analysis functions! 6 

/ Ht)$(t -n)dt = j 0(r)5(f - „) dt - 0(„) = $(„). (6J8) 

UK see in another way that^jn) = d(„), because the values of* at the integers come from 
Ijllgl eigenvector of (J, 2)2F*. When the filter is halfband, the center column of the matrix 
Mm* VeCt ° r * ™ S is 3,1 ei 8envector with A. = 1. Assuming a stable basis, there are no 
^feenvectors for X = 1 by Condition E. So 0(„) agrees with «(„) 
^^mtepiolatingproperty is highly useful in several applications. But the analysis filter 
Wpmm' ~ land *<? =/(0 is generally not acceptable. Therefore we now lift H. That 
Wgm a neWpair & • F ^ which ™nts extremely promising. 



<::■ 

I 



§ 



4 
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Biorthogonal Filters with Binary Coefficients 

A binary coefficient or dyadic coefficient is an integer divided by a power of 2. The maxflat 
halfband filters D 2p (z) all have binary coefficients. This is clear from the Daubechies formula 
(5 75) where the binomial coefficients are integers. Multiplication by a binary number can be 
executed entirely by shifts and adds. Roundoff error is eliminated. And on some architectures, 
the filter needs less time and less space. 

We are therefore highly interested in binary filters. Most factorizations of D 2p — this has 
been our route to orthogonal and biorthogonal filters-are not binary. There are zeros at irra- 
tional points like z = 2 - 73. But we can certainly move zeros at z = -1 between analysis and 
synthesis. This operation we call balancing. 

Moving from F(z) to H(z) maintains binary coefficients and symmetry. 

h neu) (n) = \[h M (n) + h eU (n - 1)] and /^(n) = \\f oli {n) -f new {n - 1)]. (6.79) 

Note/„ eio at the end. We are dividing F(z) by (^^The product F Beiu ff neu , equals F ol *Ho ld , 
so biorthogonality is preserved. The scaling function <f> new {t) isJ M (t) convolved with the box 
function. Therefore it has exactly one more derivative than <t> old «) (Section 7.3). Similarly 
S neu> (,) from Fntv, has one less derivative than <p M (t) from F oli . In a filter bank we avoid the de- 
structive factors 75, by putting both of them into synthesis. Our convention below is (1) = 1 

^C^exmple hi = [1] and/7 = [-1 0 9 16 9 0 - 1]/16 is binary. This 1/7 filter ban! 
has denominator 16. Balancing will produce 2/6 and 3/5. still binary and symmetric: 

hi = [1 l]/2 and fS = [-1 1 8 8 1 - l]/8 with 1/3 zeros at it 
&3 = [1 2 l]/4 and fS = [-1 2 6 2 - l]/4 with 2/2 zeros at jr. 
These are very effective short filters. They are probably the best-after reversing the secoru 
pair to 5/3. In the experiments of Chapters 10 and 11 they are comparable and quite effect™ 
As factors of the maxflat D 6 filter, we have seen them before. An interesting feature of /5 i 
that its scaling function 0(r) is infinite at all binary points! Section 7.3 confirms that this 00 
nevertheless has finite energy (and even 0.44 derivatives in the energy sense). By removing tw 
zeros at z = - 1 from /7, we have stolen its smoothness. Enough is left to make it good amon 
short filters (but moved to the analysis side). 

For serious compression we need more zeros— which means longer filters. Our previoi 
route was to factor a long maxflat filter. When one factor is F(z) = the pair is st 

binary and symmetric. The seating function 0(r) for this factor is a spline. It has maximu 
smoothness for its length (p intervals) coming from a maximum number of zeros at n (p zeros 
These are outstanding filters, when we keep enough smoothness in the other factor. Taking o 
three zeros from D 6 would leave a filter [-1 3 3 -l]/2 which has one zero but negative smoot 
ness-too risky to iterate. Taking three zeros from D s is allowed. Now, instead of factoring 
long filter, we will lift a short filter. 

Lifting will not maximize the number of zeros at n (although we like those zeros). Our ft 
lifting will go from W to 9/7, and we choose S(z 2 ) in (6.74) to give two zeros at jr. Here is t 
result of lifting ff = 1 (all filters are symmetric): 

ft9 = [l 0 -8 16 46 ...]/64 and /7 = [-1 0 9 16 •]/16 (2/4 zeros). 
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These are binary filters. Balancing the zeros by (6.79) yields 10/6 from 9/7: 

ftlO = [l 1 -8 8 62 ...]/128 and/6 = [~l 18 8 1 - 1]/8 (3/3zeros). 

This 10/6 pair gives better compression as 6/10— reversing analysis and synthesis. Sodoes5/ll 
from 11/5, after another balancing (or unbalancing!) step: 

All = [1 2 - 7 0 70 124 - ■ • ]/256 and /5 above (4/2 zeros). 

Figure 6.9 shows scaling functions $(t) and <P(t) for analysis and synthesis. You can see how 
the zeros affect the smoothness. You cannot easily see which pair is best in compression— that 
depends on the image. 



Analysis Scaling Function 




Analysis Scaling Function 



-0.2 




Synthesis Scaling Function 




Synthesis Scaling Function 




Figure 63: Scaling functions for A9//7 and h6/fl0. 

Note 1 The reader might be interested in the construction of these new (1995) binary filters. 
The first author found the 9/7 pair in September, but not by lifting. With /7 fixed, he solved 
the halfband equation (| 2)(f 7 * A9) = 6 for the symmetric filter h9 with two zeros at n. 
(The first zero determines the middle coefficient from the others; the second zero is automatic 
by symmetry.) When reporting this result for the Wavelet Digest, he learned that Wim Sweldens 
had created a whole family of binary symmetric filters earlier in 1995 by lifting. We propose to 
call them "binlets". Here are the next filters M3//7 and 13A//11. All signs indicate that A13 ffj 
is the right choice: 



M3=[-l 0 18 - 16 -63 144 348 -.]/512 with /7 (4/4 zeros) 
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l 3h = [-3 0 22 0 - 125 256 724 256 - 125 0 22 0 - 3]/1024 with 

/11 = Z> 6 = [3 0 -25 0 150 256 150 0 -25 0 3]/256 (2/6zeros). 

Extra length gives more zeros and higher compression, up to a point. Then ringing destroys the 
image quality. See Sections 10.1 and 11.2 for the artifacts that plague image compression. The 
boats in Figure 7.4 offer a visual comparison. 

Note 2 [Majani2] emphasizes the importance of a reversible integer implementation. Integer 
inputs are reconstructed exactly. Orthogonal transforms seem not to be reversible (except Haar 
for Af = 2 channels and Hadamard for certain Af > 2). The 2/6 biorthogonal transform with 
hi = [1 1] is reversible and very useful. Lowpass components yx come first: 

yv = *2i + JC2i+i and then y^+i = x&'+i - [yn /2J + Myv-i - yx+i)/ *6J • 

The inverse also has "even = even + /(odd)" and "odd = odd + g(even)": 

X2i = y* - *2i+i and then xy+\ = ya+i + L)W 2 J - l(yv-2 - y2i+2>/16J. 

Majani has shown that the new binary 9/7 and 13/7 transforms have reversible forms (lossless 
in integers). The normalized DCT is not reversible for integer data. 

Note 3 The maxflat Daubechies filters with 4/7-1 coefficients and 2p zeros are also known as 
Deslauriers-Dubuc filters [DesDub, CDM]. They interpolate because they are halfband. They 
leave the values x(n) unchanged and produce midpoint values x(n + 5). Section 5.5 confirmed 
that all polynomials of degree less than 2/? are interpolated exactly. Recursive subdivision start- 
ing fromx(n) = 6(n) converges to the scaling function 4>(t) by the cascade algorithm! 

Problem Set 6.5 

1. Double-shift orthogonality of lowpass filters is 2 £ *o(*)/o(* + 2/i) - 6(n). Show that in 
frequency this becomes 

H 0 (o))F Q (ci>) + H Q (co + 7r)F 0 (a> + n) = 1 . 

Write the same equation in the z-domain. 

2. Problem 1 involves a row and column of the modulation matrices F m and H m : 

F m {z)H m {z)-[ Fd _ z) Fi( _ z) HM Rx[rz) J. 

Which row-column multiplications correspond to which equations (6.56H6.58)? 

3. Suppose Ho(o))F 0 (co) + H Q (a> + ji)F q (co + n) = 1 as required. By alternating flip 

HiW^e-^Foico + n) and Fi(a>) = -e~ ia H 0 (a> + n). 

Show that the other entries of F m (z)H m (z) = / are then correct 

4. What wavelets come from the biorthogonal filters with H 0 = 1, F 0 = \z + 1 + {z'K #i = 
|z - 1 + \z~\ F\ = -1? Recognize the delta and hat 

0(0 = 2^(20 and 4>(t) = {<f>(2t + 1) + <p(2t) + \4>Q>t - 1). 
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Then construct wavelete ftom S(r) = -10(2/ + 1) + 0 (2 ,) - ifot - i) and w(t) = 
2^(2/ - 1). Check the biorthogonality conditions 

j<Kt)W-k)dt = Jw(t)w(t-k)dt = 6(k) and 
J 4>(t)w{t-k)dt = Jfct)w(t-k)dt = 0. 

have the form (6.75). The difference D = H* — H must satisfy 0(z)F(z) = D(- z )F(- z -> 
from (6J4) Substitute Z)( z ) = a + fa- + cz -> + dz -3 + . / to it 

6. Lifting (// F)to(/^^^doesnotchangethescalmgfunction0(t). Show that the new wave- 
let is w'(t) = w(t) - £ j(Jfc)0(f - Jt). 

7. The/asr*vavefe/fran^iTnissubbandmteringofmeiimerproductsa« Z.MW The 
highpass channel produces b n = </(,), S yt (/)) by (6.70): 

"^I^W-^H+U and *y* = X)*i( , - 2 *K+u- 

| Supposed is lifted to H*. Show that the lifted a%are a% = a }i + Z*{1 - k)b Jt . The fast 

m: • tnmsfbnn unhfts by -*<l - *) fa the same way. Then it inverts the <», »,) transform 

§ as usual by a/+w = £/(* - 20a* + £/",(*- 21) V 



lift 



8. From the input *(«) = a , t compute even samples a* = a,. a and odd = - 
(a* + ao.t+i). Then lift to «& = o« + ft,.*., + b^/A. Combine to recognize the 5/3 filter 
bank, computed more efficiently and in place —no auxiliary memory. 

9. Which filter h gives linear interpolation at each step of recursive subdivision? What is 0(0? 
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Wavelets and Dilation Equations 

Gilbert Strang Siam Review 31 (1989) a3 _ m 

Abstract 

Wavelets are new families of basis functions that yield the representation /(*) = 
£ 6y W(2 J x - ft). Their construction begins with the solution <f>(x) to a dilation 
equation with coefficients c t . Then W comes from <p, and the basis comes by 
translation and dilation of W. It is shown in Part 1 how conditions on the c t lead 
to approximation properties and orthogonality properties of the wavelets. Part 2 
describes the recursive algorithms (also based on the c t ) that decompose and re- 
construct /. The object of wavelets is to localize as far as possible in both time 
and frequency, with efficient algorithms. 



Wavelets are based on translation (W(x) -+ W(x + 1)) and above all on dilation (W(x) -> 
W(2x)). It is remarkable how long it has taken for "dilation equation?' to be mentioned beside 
differential equations and difference equations. True, they are hardly in the same league. But 
ideas about wavelets are coming fast The mathematics is attractive and several important ap- 
plications seem to fit — I hope this survey will be helpful. You should know that its author is 
neither an expert nor an evangelist 

The goal is a new way to represent functions— especially functions that are local in time 
and frequency (or space and wave number). Compare with Fourier series. Sines and cosines are 
perfectly local in frequency, but global in x or t . A short pulse has slowly decaying coefficients 
that are hard to measure. To reconstruct the pulse, a Fourier series depends heavily on cancella- 
tion. The whole of Fourier analysis, relating properties of functions to properties of coefficients, 
is made difficult (some say interesting) by the nonlocal support of sin x. 

In achieving local support we lose the greatest property of the basis {e inx }. With respect to 
a wavelet basis the differentiation operator is not diagonal. Wavelets are not eigenfunctions of 
d/dx, and frequencies are mixed up. The uncertainty principle imposes limits on what is possible 
m x and £ together. The commutator 0/3x)(3/3£) - (3/3f )(3/3;c) is a multiple of the identity 
(since (d/8x)(xu) - x(du/dx) = u), so we cannot diagonalize both operators. But a good 
"microlocalization" leaves 3/3* nearly diagonal, and at the same time nearly diagonalizes 3/3£ 
(which is multiplication by x). To connect dilation with multiplication by x, differentiate f{cx) 
with respect to c at c = 1. 

The second important property of {e inx } is orthogonality. That can be saved. Wavelets can 
be made orthogonal to their own dilations and translations. Then / W(x) W(2'x -k)dx =0 for 
all integers j and k. The wavelet basis has two indices, in which k is translation and j is dilation 
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or / in which b Jk carries information about / near £ = V and r - o-Jt -n, iUlul *° n 
the detail at the scaling level h = 2"/. ~ km ™ e Sum on * 18 

Orthogonality is not easy to achieve with local support. Truncated at 7m w ->„ 

to r^ ^ f UPP ° rt ty briDgin8 <** g^^ytozero.butitisnotfuUyltisr. 

to^ The price of orthogonality with compact support is irregular basis functions. We 1 v^h 
these wavelets by doing ail computations recursively (this subject is recursion heaven TlZt 
is important to recognize that orthogonality and even linear M^l^V)^ ^esselS 
in the representation of functions. Wavelets need not be orthogonal 
tn ^ bnef , m ^ ucti ° ncann ot do justice to the applications. Nor can we attempt a proner his 
tory-tf would be mostly in French. The idea of wavelets grew out of seisn^c anZ^; 
development been led by Yves Meyer, whose book will Scribe a neTZSfntoc 
andy^ (connectmg to work of Calderon, Grossmann, Morlet, Coifman, We5ta maT 0 m 
ers) The interest m wavelets is both puce and applied-like the interest in splines 

ri™ f! IJO* ^ * e P 10 !**** of wavelets - approximation through Condi 

tion A and orthogonality through Condition O. Since we never si wavelets a^ ZL 
recursively), their properties have to be discovered indirectly. W?ZSv neS 
erties in order to have any idea what the algorithms are producing lhTp J?L^ *T 

^^uatiTns^ 



1. Dilation Equations: Construction of 0 

The basic dilation equation is a two-scale difference equation: 



(A.1) 



We look for a solution normalized bv V<A/fr — i Tk»««.« 

comes rrommultiplyingby2an7m te gTll: " ^ e fot »e coefficients c t 
2 f<t> dx = T l Ckj4>(2x-k)d(2x-k) yields J>=2. 

e^Xset clt^^ * ^ = * * * not — ** * striking 

The delta function 0 = 5 satisfies *(*) = 25(2*). 
Th^diiation of. is unfamiliar (but somehow very pleasing). For other c's, spline functions 



<Kx) 



Box function: 



4>(2x) 



<H2x - 1) 



c 0 = 1 
c, = 1 
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Hat function: 




We now outline three constructions of the "scaling function" 0. Those constructions display 
very clearly the mathematics of dilation. Then we turn to wavelets, their properties and their 
purpose. A wavelet W(x) is a second combination (involving the same recursion coefficients 
c t ) of the translates <f>(2x - k). 



Collection 1. Iterate = £ c t ^-.,(2x - k) with the box function as When 
c 0 = 2 the boxes get taller and thinner, approximating the delta function. For c 0 = c, = 1 the 
box is invariant: <pj = For ±, 1, \ the hat function appears asj-^ooandJ-^S* 1 
yields the cubic B-spline. An example that will be important (an inspiration of Daubed^— we 
propose the notation D 4 ) has coefficients i (l + VJ), I (3 + VS). I (3 - and i (l - VS): 




TTus scaling function i> 4 leads to orthogonal wavelets. It is not as smooth as it looks Note 
that the Weierstrass nowhere differentiable function, which is£fc» cos(3"jc), involves dUation 
by 3. So does de Rham's function, which has c* = §, 5, 1, |, | adding to 3. Resnikoff has found 
a connection between Weierstrass functions and wavelets. 



Construction 2. The second construction takes the Fourier transform of (1): 
0(f) = £c t J<P(2x-k)e l(x dx 

The symbol P($) = i £ Ct e'*« is the crucial function in this theory. Note that P(0) = 1 Now 
repeat (2) at £/2, £/4, ... and recall £(0) = f<j>dx = 1: 

^ ) = [n P (^)]^) approaches f[p{ij). (A.3) 
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For c 0 = 2 we find P = 1 and £■ 1, the transform nf «t. *u - . 

products of the P's are geometric series: de,ta ""ebon. For c 0 = n = 1 the 

- id + ««' 2 ) (i + = Lzf!! 

7 4(l.-eW<)' 

As N — ► oo this approaches the infinite product (\ — i$\ x 
form of the box function. The hat function comes frn * ^ ™ S is & dx> ^ trans " 

(Multiplication of P's is i times convolution which b y < 3 ) 3180 s <H*res 

squaring again. * or e's). The cubic B-spIine comes from 

Construction 3. 

at the integers x 

quarter-integers, i ^i ULS 

wavelet calculations use recursion. * ** " IU * * "™ */2 y . This is fast to program. All good 

The values of <p at the integers come from an ei<**» 
ficients, set x = 1 and x = 2 in the dilation eoJE n ^ tor - Wilh *e four Daubechies coef- 
0 < x < 3: ^ Uabon (D and use the fact that 0 = 0 unless 

= H3 + V3) 0(1) + Kl + ^ (2) 

- .Hi -^*<l) + I(3-V3)0(2) (A ' 4) 
This is <f> = with matrix entries L i} = c Z ' 

ference equation. The eigenvalues are 1 and ± Th . mpare Q-j for an ordinary dif- 
#0) = * (l + V3) and ,<2) = J (l - Vs\ which *" * = 1 haS 

eigenvalue X = I means that the recursion can be d^ ^ " ^ ^ ° f ^ ^ 
leads similarly to 0'(1) and 0'(2). In some weak sen^T* 1 *** *' (x) = EcttyQx - k) 
For the hat function, the recursion matrix (see below* ' ^ = ^ has a "dilational derivative." 
the eigenvalues are 1, \, 1,1 ' agam has A = L 5- For the cubic spline 

To repeat for emphasis: From fl>(l) and 0 (2) the re^ 

In these constructions the properties of P(£) = r ?^ l °n gives everything. 

potheses are in flux, and infinitely many c t can be all 5 ^ m decisive - 1116 precise hy- 
gether the theory of dilation equations, before we on ° ° ne basic P^t^y will bring to- 

«° 00 to wavelets. 

1.1. Dilation equations: Fundamental theorem, Th 

approximation, by splines or finite elements, depends if accuracv of piecewise polynomial 
degree p ± 1 can the polynomials 1, x, x 2 , . . . , x p~i L? n e a* 8 ™* <» this question: To what 
ing functions? When the polynomials are "in the spac^' I h Pr ° dUCed exactlv bv * e approximat- 
In our case, the approximating functions are and'-, ap P roximation err °r is of order h p . 
proximation, and finite elements have the narrowest translates - Splines are the best at ap- 
we require orthogonality. There is already a theory f Upport ~ but oo* are weeded out when 
P with the properties of The link is the Poisson s^n^ 0 *"**** 011 by trans,ates - " connects 
equation, that throws new questions into the meory--? • formula - 4> solves a dilation 
questions have the same answers. ls extren »ely satisfying that these new 

For approximation with accuracy h p , the Fourier tr t ~ 
at all points £ = Inn (except at £ = 0 where 6 = ,r^ 0nn * must ^ve zeros of order p 
condition on the symbol P. According to (3), the transf ^ tiCe h ° W easUy Aat converts «° a 
At f = 2n the first factor is P{n). At £ = 4jt the seco 0 ^^ ™ 1,16 infinite P roduct of ^ 
fira factor is P(3jt), which by periodicity is the same becomes P ( 3r >- At ? = &r the 

of 0: 38 ^ Or). The zeros of P produce zeros 
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Condition A. The symbol = has a zero of order p at f = n . Equivalently, the 

coefficients ct satisfy the sum rules that yield P (m) (^) = 0: 

^(-l^c^O, m = 0,1,..., p-l. (A . 5) 

The box function has P = \ (l + e*« ) and /> = 1. The hat function has p = 2 and so does 
D4. The cubic spline has p = 4. 

A zero at £ = ;r/2 (instead of ji) would also produce the desired zeros in the product £ 
Thus Condition A is not strictly necessary in what follows. Choosing c 0 = 1 and c 2 = 1 and 
p ~ 2 0 + c2,{ ) stretches out the box function— it becomes <f> = \ on the double interval 

0 < x < 2. But />(*/2) = 0 produces instability and linear dependence— the alternating 
sum of stretched boxes is E(-1)*0(* - k) = 0. With the added requirement of stability, the 
condition is exactly right 

The fundamental theorem states the consequences of Condition A: 

1. The polynomials 1, x, . . ., jc*" 1 are linear combinations of the translates 4>(x - k). 

2. Smooth functions can be approximated with error 0(A") by combinations at every scale 
h = 2~ J : 

\f ~ X>*(2'* < C2"" fl/WJ for suitable a k . 

k 

3. The first p moments of the wavelet W(x) (see below) are zero: 

Jx m W(x)dx = 0 form = 0,...,p-l. 

4. The wavelet coefficients of a smooth function decay like | / f(x) W (2 j x) dx\ < C2"^>. 

5. The recursion matrix M N determining <f> at the integers has eigenvalues 1 , \ , . . . , (I) p ~ ! . 

1 and 2 come from approximation theory. The combination of 0's at scale j is also a combi- 
nation £ b jk W (2'x - k) down to scale j. 3 and 4 are easy once wavelets are defined. Mallat 
gives a sharp result, with properly stated requirements on the smoothness and decay of <p: The 
HP norm of / is equivalent to the corresponding norm of its coefficients b jk . Wavelets lead to 
unconditional bases, suitable for a wide range of function spaces. 

It is 5 that makes 4>(x) smoother as p increases and also makes the constructions successful 
The smoothness is weaker than 0 e C*-\ but it is striking that "dilational derivatives" come 
at the same time as higher degrees of approximation. What remains to be studied is orthogonal- 
ity — which imposes an entirely different condition on the c* . 

Remark i. Suppose the basic recursion has coefficients c 0f . . . , c N . Then <f> is zero outside 
the intervals N]. With continuity it follows that 0(0) = Oand0(tf) = 0. Those were assumed 
in (4) when we determined <p = D 4 at the integers. For the box function with N = 1, 0(0) 
and 0(W) cannot both be dropped. Our recursion matrix will be (Af^) 0 - = cy-j with i, j = 
0, . . . , N - 1. For the box function M x - [1] has eigenvalue X = 1, as expected in S above. 

The spectrum of the infinite matrix M (allowing all i, j) is an attractive problem in operator 
theory. Notice that M is convolution followed by decimation— multiplication by the matrix c w 
followed by projection onto even-numbered coordinates. By contrast with the usual Toeplitz 
case, eigenfunctions can have compact support! Homogeneous difference equations with zero 
boundary conditions lead to 0 = 0, but not so for dilation equations. 
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Remark!. The minimum requirement is p = 1. Then P(it) = 0, which means that Ten = 
l^c u +i. Smce 2>* = 2, thecolumns of M add to 1: 



Co 

C 2 Cl c 0 

*3 C 2 Ci 



steps of 2 down columns 
steps of 1 across rows 
here N = 4 



(1. 1, 1, 1) is a left eigenvector with A. = 1. The right eigenvector yields the values 0(0) 
<p(N - 1) at the integers. The recursion determines <f> at all dyadic points. Values at other Lints 



are never used. 

1.2. Wavelets and orthogonality. Finally we define a wavelet It comes from the scaling 
function 4> by taking "differences": 6 

W(x) = £(-l)V**(2* - *). (A6) 

We write W in place of the usual f, to distinguish more clearly from <fi. Notice 2x on the right 
and especially (-1)*. Examples show the effect of alternating signs: 

W 2 (x) 




Haar wavelet from box function "Wavelet" from hat function 

W 2 (x) = <f>{2x) -4>(2x-l) W = 4>(2x) - I^(2x - 1) - i<f>(2x + 1) 




1-5 2 2JS 3 

Wi(x) from (j> — D^ Orthogonal wavelet 
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The wavelet from the hat function does not belong here. It is not orthogonal to W(x + 1) 
The point is that the other two do belong. The Haar function is orthogonal to its own translations 
and dilations. Historically it was the original wavelet (but with p = 1 and poor approximation) 
The orthogonal wavelet W 4 has p = 2 and second-order approximation. 

Without formulas for D 4 and W 4 , how is the orthogonality of their translates known? We 
need a test that applies to the recursion coefficients c t , or to the symbol />(£) = I £ Cte «'*f 

Condition O. 

\P(l)\ 2 + \P^ + n)\^l or J^e^mTg^. 

With this condition, the infinite matrix L"L in Part 2 is an orthogonal projection To see now the 
role of Condition O, suppose the functions 0»(2* - *) are orthogonal. Then so are the translates 
of fa (x) = £ c k fa{2x - fc): 

f<Pi(x)<h(x-m)dx = /(E^(2*-*))(X)^(2*-2m-/))^ 

= £c*<*-2m f$(2x)dx = 0 for m ^0. (A.7) 

Construction 1 creates <p by iteration from the box function, which is orthogonal to its translates 
Therefore (as Daubechies observed) so is <f>. 

The wavelet W(x) in (6) is also orthogonal to <f>(x-m). This is simple but neat, not involving 
Condition O. The sum in (7) changes to 

y^(-l)*ci- t ct_ 2w . which is identically zero! (A.8) 

Just replace* by l-n+2m. This identity is HL* = 0inPart2. Then (6) makes W(x) orthogonal 
to W(2x - m). The orthogonality of W(x) and W(x - m) comes back to Condition O. 

The goal in constructing wavelets is to satisfy Conditions A and O. The basic family W 2 , W 4 
W 6 ,... was discovered by Daubechies, following Haar's W 2 . The accuracies are p = 1, 2, 3 
and there are 2, 4. 6, . . . nonzero coefficients c k . The smoothness also increases with' p— but 
only by about 2 derivative each time. D 4 and W 4 are Holder continuous with exponent .550 • • 4Wf 
In Galerkin's method for solving differential equations, it is natural for these wavelets to be tbf v ;i <■ 
trial functions— broader support than splines, nonsymmetric but orthogonal, multigrid built in; : &k 
all computations based on recursion, difficulty to be expected at boundaries. The first experi- 
ments by Glowinski, Lawton, and Ravachol are particularly interesting for Burgos' equation^ |J|j 

2. Algorithms for wavelet expansions 

Now comes a change of direction. Instead of discussing the properties of wavelets, we described j 
algonthms. The main question is how to decompose a signal into its wavelet coefficients, and " 1 
how to reconstruct the signal from the coefficients. There is a "tree algorithm" or "pyramid al- v ^^ 
™«*w»*u„* — ^ . • ■ ■ «- - - rJ :r -^M 



gorithm" that makes these steps simple and fast It does for the discrete wavelet transform wharf 

tfl<» Pact Pntin'or TrmncfAm /T3L."1*\ J_ r *1 l? . «. - . ^ Vj 



— — r~ ****** « uwo iui uic uiauicic wavciei uaiiMunii wuat 

the Fast Fourier Transform (FFT) does for the discrete Fourier transform. The algorithm is fully 
recursive. " J ' "^Splgf 

The user chooses a specific wavelet We begin with the simplest choice, based on the bral|M 
function. It satisfies the orthogonality property (Condition O), so all pieces of the decompos|i^ 
tion are orthogonal. The approximation property (Condition A which preserves polynomiabftS 
determines how quickly the coefficients decay— for efficiency we want to stop the decbm^l^ 
sinon early. In that respect the box function is poor. Efficiency is the reason for working wiu||l|| 
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-and simplicity is thereasonfor starting with W t . ThisisHaar's 



higher wavelets W 4 , W 6 , W 8 , . . 

wavelet [1 -1]. _ u 

/ = /* + /«» + . 

The "detail" /t/> i s a combination of 2' wavelet* a * «.„i. i- / j ^ . 

Action For a numerical exam^ taT / I 2 ^ L fi^ ^ f ^ ° f ^ 

terms, here with coefficients -4 and ft„ - ^ finestdetail /<»> „ the sum of two 



/ = 



9 
1 

2 
0 



= 3 



1 
1 
1 
1 



+ 2 



1 
1 

-1 
-1 



+ 4 



1 

-1 
0 
0 



+ 1 



0 
0 

1 

-1 



(A.9) 
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: 9 " 




" 5 ' 
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2 




1 




0 




1 



(A.10) 



This is followed by a permutation, in which high frequencies go to the bottom: 



1 




" 5 " 




" 5 " 


1 




1 




i 


1 




1 




4 


1 _ 




4 




1 



The next step is another butterfly, on low frequencies only: 



(A.11) 
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2 



1 

? 

~2 



5 
1 
4 
1 



3 
2 
4 

L 1 J 



(A. 12) 



The result is the set of wavelet coefficients! 1 9 4 i tk 



l 

4 
I 
4 
1 

2 

0 



1 
4 
1 
4 
1 
2 

0 



1 
4 
_! 
4 

0 

I 

2 



1 
4 
1 
4 

0 

1 

"2 



has D~ l = R = 



1 
1 
1 
1 



1 
1 
-1 
-1 



1 0 

-1 0 

0 1 

0 -1 



The coefficients 3, 2, 4, 1 enter the vector b - (h u u u x ^ 
in (9) is / = * bm ^ coefficients »7^7> fc? ^ CXpansion 

/c / - £>/. This product £>/ was computed 
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recursively, from two butterfly matrices with a permutation between. In general there will be J 
matrices with permutations between. 

The reconstruction is also recursive. It inverts (12) then (11) then (10). The global matrix R 
is the product of these local inverse matrices. 

Notice that the operation count is proportional to n. It is best possible (the FFT count is 
n log 2 n). There are only n - 1 individual 2-by-2 matrix multiplications, since high frequency 
coefficients (here 4 and 1) are setded and not reused. The Walsh functions give a different piece- 
wise constant representation, in which the last two basis vectors are (1, -1, 1, -1) and (1, -1, 
-1,1). In that case 4 and 1 enter another butterfly to produce the Walsh coefficients § and y 
The Walsh basis is global. The wavelet basis is local, but scaled — its support has width O (2~ ) 
at the finest scale and O (1) at the coarsest scale. 

Notice also the normalizing factors \ in decomposition (and l's in reconstruction). The al- 
ternative is to introduce 1/V2 for both. This has the advantage of normalizing the wavelets 
Wjt = 2 i/2 W(2 J x - k) at every scale. The whole basis is orthonormal (when || W|| = 1). 
In the discrete case R and D become orthogonal matrices: 

has R = D" 1 = transpose of D. 



Based on the Haar example, we now start on Mallat's beautiful tree algorithm for wavelets. 
The simple averagefrom [\ \] is replaced by a discrete filter based on 0. The difference [\ -5] 
is replaced by a filter based on W. The filters use the same recursion coefficients c k that led to 
<f> and W in the first place. 

Decomposition. The given n-vector / is on the finest scale h = 2~ J . The fine-to-coarse fil- 
ter (the "restriction operator" in multigrid language, the lowpass filter in signal processing lan- 
guage) is L. It produces a vector with half as many entries: 

(iA-jE^-i/i. * — 1 2* (A ' 13) 

In the Haar example with c 0 = c x = 1, the entries of Lf are \{f\ + fi) and |(/ 3 + /*). 
The recursion continues to coarser scales, and after / steps it reaches a single number — the 
coefficient b+ in /* at the coarsest scale h = 1. Here b# = \ (f\ + h + h + U)- 

The dual to L is the coarse-to-fine map L* (the "interpolation operator" in multigrid lan- 
guage). Notice the change of index and the disappearance of \: 

(L*g) j =Y f C 2i - jgi , ; = l,../,n. (A.14) 

In the Haar example L*Lf has entries £(/, + / 2 ), J(/i - / 2 ), \(h + f*)> \(h ~ /*)• li is 
the projection of / onto the subspace that is piecewise constant at scale 2ft. It gives a blurred 
picture, with details lost 

The decomposition picks out these details, orthogonal to the average. The projection onto 
the wavelet subspace is the high frequency component: 

fV~» = f-L*Lf. (A.15) 

This repeats at every stage. There is an "average" or "blurred picture" a u " X) - La u \ starting 
from a U) = /. The detail lost in that average is the component of / at that stage: 

f U-i) = (/ _ L * L ) a W = a U) - L*a u - l \ (A.16) 
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^taST Statement ° f Ae dCCOm P° Siti0n d g° rithm - We will see how Condition 0 simplifies 
Reconstruction. To produce /from its details /«, run the xecmsion (16) inreverse: 

«°- , = /°- ,, + raO-). (A17) 

Ita starts from the coarsest detail /«» and the totally blurred picture a™ = ft ft returns to 
Review the Haar example first. The lowpass filter gave a") from / = a< 2 > : 



L ci Co \ 



fx 

h 
h 
n 



9 
1 
2 
0 



-[ ; ] 



ffic^ffiTnC = (5 ' 5,1, ! ) - At * en ^tevel the low-pass filter leaves 3, the co- 
Seifp^uces ^^^^^onalfilter-thehighpass filter H. In the Haar 



"U co -c,J 



fx 

h 

h 
L f* J 



9 
1 
2 
0 



Z^f^f ^ ^ I T 5Sent 1,16 detaa 7 - (4 ' ~ 4 ' »• which when a« 
. blurred to ^ At the next level is applied to «« That produces 1(5) - 1(1) = 2 This 

the coefficient b 0l , representing the detail (2, 2, -2, -2) lost when J> is bluned to > We 

ow put these pieces together into Mallat's pyramid algorithm- ' We 



womporftion. Initialize a' = /. For j = 7, . . . , l compute 

a J ~ 1 =La i and b*- 1 =HaK 
econstruction. Start with *° and 6°. . . . . For j m j y 

le full decomposition is represented by a tree of filters: 



(A.18) 



(A. 19) 



b J ~ x 



a'~ 2 
b'~ 2 



b° 



ie reconstruction goes from the branches of the tree back to the root 



a" 
b° 



b l 
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The next step is to identify these filter matrices L and H for examples other than "box and 
Haar." 

Note. The filter matrices L and H have half as many rows as columns. By dropping the 
parentheses around /, we distinguish the vector a j with only 2 s components from the vector 
a U) with the full w J = n components. The vector a j contains the expansion coefficients of a U) 
with respect to the translates <l>(2 j x - k). See the example above and the multiresolution below! 



2.1. The filter matrices L and H. The matrix L is known from the first part of the paper. 
Its entries L t j = ca-y are the recursion coefficients for the scaling function. Rows 1, 2 and 
columns —1, 0, 1, 2 are displayed with N = 3: 

L = - T C3 Cl Cl C ° 1 
2 [ c 3 c 2 ci c 0 J* 

The beautiful thing is that the highpass filter (strictly speaking it is band-pass) uses the same 
coefficients. H is associated with the wavelet W just as L is associated with the scaling func- 
tion (f>. Equation (6) for W uses the same c*, but with alternating signs and reversed order. The 
wavelet filter has 

fliy=(-lV +, c i+1 _ 2f . (A.20) 
Rows 1, 2 and columns 1, 2, 3, 4 are displayed: 




The indices were chosen to match the Haar example (variants are possible). The transposed ma- 
trices, without the factor \ , represent the dual filters L* and H* . Hie important points now come 
quickly, and matrix multiplication is the best proof. 

Theorem 1. By their construction the filters are orthogonal: 

ffL*=0. (A.21) 

This multiplication is the reason behind the construction of H — alternating signs, reversed or- 
der, index shifted by one. See equation (8). 

We finally come to the reward for Condition O: £ cjfeCt +2m = 28om. The reason for that 
condition is in the reward. Remember that the box function and D 4 satisfied this requirement but 
not the hat function or the cubic spline. Condition O can be stated and understood in transform 
space, but I believe that the matrix interpretation is again the clearest 

Theorem 2. If condition O holds then 

1. LL* = / and HH* = L (A.22) 

2. L*L and H*H are mutually orthogonal projections with 

L*L + H*H = /. (A.23) 

Remember that L and H map into subspaces half as large as the original. L* and H * map back. 
The identity operators in (22) are on the half-sized subspaces. 

The proof of (22) is by direct matrix manipulation. Condition O gives the result. Then it 
follows that L*LL*L — L*L t so L*L is a projection — and similarly for H*H. The property 
HL* = 0 in (21) yields H(L*L + H*H) = H. The transpose LH* = 0 yields L(L*L + 
H*H) = L. The operator in (23) is the identity on both orthogonal components — the ranges 
of L and H — so it is the identity. We have an orthogonal decomposition by "quadrature mirror 
filters*' L and H at every step. 
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2.2. Multiresolution of L 2 . The last paragraphs changed quietly from functions to vec- 
tors. That was for the sake of algorithms, which use values of <f> and W at dyadic points k/2 J . 
The Haar example began with / at equally spaced points on (0, 1]. But the filter matrices really 
apply to discrete values along the whole line— they are infinite matrices. More than that, the 
decomposition / = £ fW is just as valuable for functions in 1? as for vectors in I 2 . 

This multiresolution yields the details of f at all scalings 2~ J . On the whole line we take ; = 
0, ±1, ±2, . . The decomposition develops an idea that was already present in approximation 
theory — to put frequencies together in "octaves." (Besov spaces combine frequencies 2 s < 
£ < 2 j+1 . It seems that the ear also receives frequencies on a logarithmic scale.) For functional 
analysis the starting point is the subspace Sj spanned by the translates </>(2 j x - k). If a function 
g(x) is in Sj , then g (2x) is in Sy+i . The dilation equation writes <j> (x) as a combination of <f> (2x - 
*), which assures that 5 0 C Si . At all scales we have 

■••SLiCSbCSi C5 2 --- with U Sj dense in L 2 and CiSj = {0}. 

Nov/ turn to the wavelet subspaceWj . It is spanned by the translates W(2 J 'x-k). It is invari- 
ant under translation by multiples of 2~K If g(x) is in Wj then g(2x) is in Wj+ X . The construc- 
tion W(x) = E(-l)*Xi-*0(2* -*) puts Wand itstranslates into S u and makes them orthogo- 
nal to 5o. In fact, W 0 and Sq are orthogonal complements in S\ . At every scale Wj ® Sj = Sj+\ . 
The spaces Sj give the "partial sums" of the differences W): 

•••0W-i©Wo©~-®W y = S, +1 and ©W0=ZA 

The multiresolution of / is a splitting into components /U> eWji 

/ = £/ 0) or / = /* + E/ 0 " ) 1 f+eSo. (A.24) 

-oo 0 

This is a very satisfying decomposition of L 2 functions, classical but with new subspaces. The 
coefficients b J in Mallafs pyramid algorithm corresponded to f U) e Wj, and a* corresponded 
toa^eSj. 

The analogue of the discrete Fourier transform was in the algorithm. The analogue of ordi- 
nary Fourier series is (24). The analogue of the Fourier integral formula is the integral wavelet 
transform. Representations of different groups give rise to different transforms. 

2.3. Applications. Image processing works with F(x, y), so it is natural to look for two- 
dimensional wavelets. The simplest construction uses the products 4>(x)<f>(y), <p(x)W(y) 9 
W(x)<j>(y), WQc)W(y). Orthogonality is clear. New constructions have been invented that are 
genuinely two-dimensional, but it is useful to start with the tensor products of "box and Haar." 
The given two-dimensional array F yieWs a two-dimensional array B of wavelet coefficients. 

For pattern recognition, a major difficulty with the wavelet transform B is the lack of transla- 
tion invariance. If the pattern is shifted by a fraction of h, its wavelet model is changed. A higher 
sampling rate is possible but expensive. Mallat studies instead the zero-crossings of the wavelet 
transform, which locate the signal edges. Now the difficulty is to make the reconstruction stable. 
In edge detection the first wavelets were Laplacians of shifted Gaussians, introduced by Gabor. 
The orthogonal wavelets of Meyer are C°° with polynomial decay, the Battle-Lemarid wavelets 
based on splines are C n with exponential decay, and the Daubechies wavelets are C n (smaller n) 
with compact support 
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Wavelets and Filter Banks 

In closing we recall the original problem-to localize in time and frequency. Geoohvsi^ 
needs to represent short high-frequency pukes. Physics needs to divide up phase space TW 
herentstates^ = e^gOc -q) giye a "Weyl-Heisenberg" frame, with some redundancy - htt 
st.ll / can be constructed from ff(f, g„)g„dp dq. Mathematics needs (or wants) an who? 
onal decomposition, better than g„ at high frequencies and with no redundancy. The answer for 

wavelets* ^ ***** Daubechies and Howard Resnikoff for introducing me to 
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